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1.  COMPUTER  PROGRAM  ABSTRACT 


PROGRAM  NAME: 
ORIGINATOR: 


CONTRACT  MONITOR: 


CONTRACT  NUMBER: 


MCARLO 

Bolt  Beranek  and  Newman  Inc. 

50  Moulton  Street 

^|rab^i^|e , ^Massachusetts  02138 

David  L.  Kleinman,  Sheldon  Baron, 
and  Jeffrey  E.  Berliner 


U.  S.  Army  Missile  Command 
Aeroballistics  Directorate 
(205)  876-1951 
Richard  E.  Dickson 


DAAH01-76-C-0194 


PROGRAM  ABSTRACT 


MCARLO  is  a computer  program  for  generating  simulated  time  histories  of 
pertinent  variables  in  a man-machine  control  task.  The  optimal  control  model 
(OCM)  forms  the  basis  for  the  Monte-Carlo  simulation  equations.  The  ensemble 
statistics  of  such  time  functions  must  agree  with  the  covariance  propagation 
results  obtained  via  more  direct  methods,  e.g.,  using  TIVAR. 


The  MCARLO  program  is  written  in  the  FORTRAN-IV-EXTENDED  computer 

programming  language , and  is  designed  for  efficient  batch  operation  on  a 
Control  Data  CDC-6600  computer.  Data  input  to  the  program  is  provided  on 

standard  punched  cards  and  output  is  generated  via  the  linepr inter . 

In  this  manual  we  give  the  modeling  formulation  and  requisite 
discretization  of  the  equations,  a description  of  the  MCARLO  subroutines,  the 
input  deck  setup  and  a sample  problem  with  solution. 
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2.  PROBLEM  FORMULATION  AND  ALGORITHMS 


The  major  aspects  of  computer  simulation  for  a man-machine  system  are 
shown  in  Figure  1.  The  modeling  issues  are  discussed  below. 

2. 1 System-Display  Dynamics 

There  are  no  modeling  restrictions  on  the  system  being  controlled,  other 

than  the  generation  of  a set  of  NY  displayed  elements  y time  t from: 

u 

u(t)  = human's  control  inputs,  NU  vector 
w(t)  = random  input  disturbances,  NW  vector 
z(t)  = "deterministic"  inputs,  NZ  vector. 

In  its  most  general  mathematical  form,  the  system/display  dynamics  might  be 
modeled  by: 


where  x(t)=  NX  system  state  vector.  The  vector  w(t)  consists  of  independent 
zero  mean  white  Gaussian  noise  inputs  with  covariance 


For  the  special  case  of  a linear  time-varying  system,  the  equations 
(1)-(2)  become: 


x(t)  = f(t,  x(t),  u(t),  w(t),  z(t));  x^  3 x(to) 
y(t)  = h(t,  x(t) , u(t)) 


(1) 


(2) 


E[w^(t)  Wi(a)l  = Wgi(t)  5(t-<j)  i=1,...,NW 


(3) 


x(t)  = Agx(t)  + Bgu(t)  + EgW(t)  + FgZ(t) 

y(t)  = C x(t)  + D u(t) 
s s 


(5) 


(4) 


The  system  parameters,  which  may  be  time  varying,  are: 


= NX  by  NX  state  matrix 
B = NX  by  NU  control  matrix 


s 


E = NX  by  NW  noise  matrix 
s 

= NX  by  NZ  bias  input  matrix 
= NY  by  NX  state  display  matrix 
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Figure  1,  Major  Aspects  of  Computer  Simulation  for  a 
Man-Machine  System 
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^3  = NY  by  NU  control  display  matrix 

A means  for  updating  the  system  parameters,  as  a function  of  time,  must  be 
included  along  with  the  system  description.  Since  the  form  of  updating  is 

highly  dependent  upon  the  system  model,  a general  updating  scheme  is  feasible 
for  only  the  highly  structured  linear  case. 

2.2  Human  Operator  Internal  Model 


In  the  OCM,  the  human  is  assumed  to  have  an  internal  characterization  of 
the  input-output  response  of  the  system.  This  ''internal  model"  is  assumed  to 
be  linear , in  state  variable  form. 


+ E„w_(t)  + F_z_(t) 


m m 


m m 


m m 


= Vm<«  " 


where 


(6) 

(7) 


^ni(t)  = internal  model  states,  NXM  vector 

= model  deterministic  inputs,  NZM  vector 

(present  formulation  requires  z^=z  so  NZM=NZ) 

= model  Gaussian  white  noise  inputs,  NWM  vector 
E[Wmi(t)  = Wgi(t)  <S(t-a);  i=1,...,NWM  (8) 

The  model  Inputs  u (t)  and  displayed  outputs  y (t)  are  assumed  to  be  the  same 
™ m 

as  the  actual  system  inputs  u(t)  and  displays  y(t)  to  avoid  numerous 
conceptual  problems.  Thus, 


u (t)  = u(t)  and  NUM  = NU 
m 

YmCt)  = y(t)  and  NYM  = NY. 

The  internal  model  parameters,  which  can  be  time-varying,  are: 

= NXM  by  NXM  model  state  matrix 
= NXM  by  NUM  model  control  matrix 
Ejjj  = NXM  by  NWM  model  noise  matrix 

= NXM  by  NZM  model  bias  input  matrix 
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= NYM  by  NXM  model  output  matrix  for  states 

D = nyM  by  NUM  model  output  matrix  for  controls 
m 

The  choice  of  model  parameter  matrices  is  somewhat  subjective.  In  the 
general  non-linear  case , these  matrices  typically  would  approximate  the 
partial  derivatives  of  f and  h in  Equations  (l)-(2),  i.e., 

A ~ 3f 

*m  I 

In  the  linear  system  case  of  Equations  (4)-(5),  the  model  typically  would 
reflect  an  appropriate  lower  order  characterization  of  the  true  system 
dynamics.  Of  course,  a not  unreasonable  choice  for  the  model  matrices  is: 


A = A , etc. 
m s ’ 

Here,  the  model  parameter  is  assumed  to  be  same  as  the  associated  system 
parameter.  This  is  a convenient  assumption  as  it  greatly  simplifies  the 
process  of  updating  model  matrices  for  time  varying  systems. 

The  internal  model  is  used  within  the  OCM  to  help  generate  a (continuous 
time)  human  operator  control  input  via: 


u(t) 


-L. 


xm(t) 
u(  t) 


+ Lq2  Vy(t) 


(9) 


The  NUM  by  (NXM+NUM)  feedback  gains 

Lc  = [TN'^Lopt  1 Tn"^]  = [Lcl  I 

are  generated  via  auxiliary  programs  that  solve  the  optimal  control  problem 
for  the  model  equations.  The  model  is  also  needed  in  the  construction  of  the 
Kalman  filter-predictor  that  generates  the  model  state  estimate  x^(t) . 

2. 3 Human  Limitations 

The.  human  generates  on  the  basis  of  the  delayed  and  noisy 

perceived  information: 


y 


Ni[yi(t-X)]  + 


- 5 - 


i=l, ...,NY 


(11) 


Report  No . 3^63 


Bolt  Beranek  and  Newman  Inc. 


where 

= the  human's  time  delay, 

^y(t)  = the  observation  or  sensor  noise  at  time  t, 
and  N^(.)  is  the  non-linear  observation  threshold: 

x-a  X > a. 

1 1 

N^(x)  =0  1x1  _<  aj^  ( 12) 

x+a  X < -a, 

1 1 


In  a simulation  model,  it  is  possible  to  implement  the  non-linear 
observations  using  Equations  (11)  and  (12),  However,  in  a man-machine 
context  we  find  it  more  convenient  to  replace  N^(.)  by  an  equivalent  gain, 
^j_.  The  random  input  describing  function: 

^ I Y ^ 

M I X \ 

is  used.  is  interpreted  as  the  probability  that  the  human  will  respond  to 
yi,  given  its  present  value  at  time  t. 

Each  sensor  noise  v^^Ct)  is  a zero-mean,  white  Gaussian  noise  with 
covariance : 


vyi(a)] 

^i(t) 

that  contains  both  an  additive  and  a ratioed  component: 


(14) 


V°i(t)  = Vy.(t)  +7Tp^.E[y.2(t-  )]  (15) 

The  quantity  f^>Q  attentional  allocation  to  the  displayed  variable  yi. 

The  are  constrained  by: 

NY 

I f. (t)  = f = constant  total  attention  (I6a) 

i=l  ^ ^ 

^i+l(t)  = fi(t)  i=1,3,...,NY-1  (16b) 

to  indicate  that  position-velocity  pairs  are  obtained  simultaneously  form  the 
display  elements. 
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The  neuro-motor  interface  portion  of  the  model  is- given  by  Equation  (9). 

The  motor  noises  v . (t)  , i=1,...,NU  are  zero-mean  white  Gaussian,  with 
ui 

covariance : 


E[Vui(t)  v^.(a)]  = V°i(t)  5(t-a)  (IT) 

that  contains  an  additive  and  a ratioed  component, 

V°i(t)  = Vyj^(t)  + TT  Var[u^(t)]  (18) 

2. 4 Discretized  Equations 

The  implementation  of  the  human  operator  simulation  on  a digital 
computer  requires  the  discretization  of  both  system  and  model  equations. 

Given  a computer  time  step  A,  the  system  must  generate  y(k)  = y(t^+kA)  = y(t) 
from  the  inputs  u(k-1),  w(k-1),  and  z(k-1)  which  are  assumed  to  be 
piece wise- const ant  over  the  previous  time  interval,  e.g., 

u(t)  = u(k-1)  t _<  to+kA  (19) 

For  the  case  in  which  the  system  is  described  by  the  linear  equations 
(4)-(5),  the  discretization  chosen  is: 


x(k+1)  = >I>gX(k)  + TgUCk)  + A[EgW(k)  + F^z(k)] 


y(k)  = C x(k)  + Dgu(k-I) 


(20) 

(21) 


vriiere  x(0)  = x(t  ) = initial  state.  The  discrete  system  matrices  $ and  F_ 


are : 


$s  = e 


Ac*  A 


r 

> s 


•0 


A3  CT  , 

e ^ Bgda 


(22) 


Discretization  of  the  linear  model  equations  (6)-(7)  is  done  in  a manner 
similar  to  the  above.  Thus 


X (k+1)  = ^ X (k)  +r  u(k)  + A[e  w^  + F z (k)] 
m mm  m m m m m 

y(k)  = C^xni(k)  + D^u(k-I) 
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where 


■ ra 


= e 


Am  A 


m 


e B^da 


(25) 


The  human  operator  model  must  generate  a control  input  u(k),  to  use  over  the 
time  interval  (t,t-fA),  via 


u(k)  - u(k-1) 
A 


-L 


u(k-1) 


+ L2Vy(k) 


(26) 


Note  that  it  is  the  control  input  itself  that  is  considered  to  be  piecewise 
constant  for  interface  with  the  system  model.  This  is  in  contrast  to  the 
covariance  propagation  approach  where  control-rate  is  assumed  piecewise 
constant  with; 


u(k) 


(27) 


The  gains  L=[L^1L2]  in  Equation  (26)  are  computed  from  either  the  gains  or 
according  to 


L 


(28a) 


or 


respectively,  where  are  the  equivalent  discrete  gains: 


1 

fA  T 1 
eAc^da 

6 

1 

- 

K 1 

■-c=Ti 

1 

~^c2 

The  discretized  observations  are: 

Ypi(k)  = N^  yj^(k-N)  + Vyj^(k) 


(28b) 


(29) 


(30) 
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where  N=  integer[T/A]  and  the  covariance  of  the  piecewise-constant  white 

noise  v . (k)  is  a -1 [Vo. (k)/f . (k) ] to  account  for  the  finite  time  step. 
yi  ^ yi  1 

Similarly,  the  covariance  of  the  motor  noise  now  becomes  vSi/A . 

2.5  Human  Operator  Model  Equations 

Equations  (23),  (24),  and  (26)  may  be  combined  into  an  augmented 
man-model  equation,  suitable  for  Monte-Carlo  simulation.  Defining  the 
augmented  state  x(k)  = [x(k) ,u(k-1 ) ] , and  input  w = [”^,v^],  we  obtain 

x(k+1)  = $x(k)  + Tu^(k)  + Ew(k)  + FZnj(k)  (31) 

y(k)  = Cx(k)  (32) 

where  u (k)  = L x (k)  is  the  "commanded"  control.  The  augmented  matrices  are 
c 1 m 

obtained  by  rewriting  Equation  (26), 

u(k)  = (I-AL2)u(k_i ) + Auc(k)  +AL2vu(k)  (26a) 

and  combining  with  Equation  (23),  yielding: 


E = A 

'r  ' 

m 

I 

; F = A 

1 

; E = A 

■o-t  L 

- 

L 1 ^ J 

(33) 


ir  (i-al„)“ 

m j m L 

; c = 

c j 

0 T (i-al„) 

1 z 

mi 
^ 1 

The  equations  (31)-(32)  are  similar  to  the  discretized  equations  that 
occur  in  the  covariance  propagation  studies  using  the  OCM.  Borrowing  heavily 
from  earlier  efforts,  it  is  easy  to  write  the  equations  for  the  Kalman 
filter-predictor  combination  that  generates  the  state  estimate  x(k)  . For 
compatibility  with  the  covariance  propagation  modeling,  we  use  the 
a posteriori  estimate 


x(k)  = x(klk)  = E[x(k) ly^(o) , ...  ,yp(k)] 


(34) 
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The  Kalman  filter  generates  the  (a  posteriori)  estimate  of  the  delayed  state, 

p(klk)  = E[x(k-N)  lYpCo) YpCk)]  (35a) 

and  the  (a  priori)  one-step  ahead  prediction 


p(k+llk)  = E[x(k+1-N)  lYp(o) YpCk)] 


(35b) 


bY  means  of  the  usual  update  and  propagate  set  of  equations.  These  are, 
respectively; 


p(k|k)  = p(klk-l)  + v(i^) 

p(k+1  Ik)  = <|)p(klk)  + 

where  V(k)  is  the  innovations,  or  residual  sequence 
V(k)  = Yp(k)  - C p(k|k-1) 
and  the  initial  condition  is 
p(0|-1)  = given  = x(0) 


(36a) 

(36b) 


(37) 


(38) 


The  filter  gain  is(^) 

\ = k|k-lC'[C  + ^Vy(k)]-1 

where  the  update -propagate  sequence  for  generating  the  Riccati 

^klk  = <I-V>  ^klk-l 

^klk  = '^^klk*'  * * FZF' 


(39) 

solution  E is: 
(iiOa) 
(iiOb) 


and 


W = diag[w  . (k)/A,vo. (k)/A] 
mi  * ui 

Z = diag[I^J"  Z^2(k)] 

are  "pseudo-noise"  covariance  matrices.  The  "correlation  time", 
^cor  = 1 sec. 

The  predictor  forms  the  estimate  x(k)  from  p(k|k)  using: 


For  simplicitY,  the  RIDF  gain  is  included  with  Vy,  rather  than  with  C. 
This  is  the  usual  practice  in  the  OCM. 
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N-1 

i(k)  = ,N  p(klk)  + I $iruc(k-i-1)  N>0  (41) 

i=0 

If  N=0,  x(k)  = p(kl k)  . 

2.6  Special  Considerations 

There  are  several  issues  in  the  simulation  of  the  above  equations  that 
remain  to  be  resolved.  Some  are  unique  to  the  man-machine  problem. 


2.6.1  Storage  of  Delayed  Quantities 

The  simulation  equation  (41)  requires  knowledge  of  the  commanded  control 

u^(k-1 ) , . . . ,u^(k-N)  . The  update  equation  (36b)  requires  u^(k-N).  Similarly, 
the  computation  of  v(k)  is  based  on  the  delayed  quantity  y(k-N).  To  meet 

these  requirements  we  retain  in  storage  , at  time  k, 


PASTUC  = [u^(k_N)  ....  Uc(k)]  (42a) 

PASTY  = [y(k-N)  y(k)]  (42b) 


2.6.2  On-line  Variance  Estimation 

The  diagonal  observation  noise  covariance  matrix  Equations  (39) 

and  (40a)  is  given  by: 

V Vy(k) 

\i(k)  = 1 i=1,...,NYM  (43) 

fi(k)  N.2(k-N) 

where  f (^)  is  the  fractional  attention  to  y . at  time  k>N , and: 
i pi  — 


^yi(J=^)  = Vyi(k)  + PyiE[ yi^(l<“N ) ] 

^lyj(k-n)  I 


N. (k-N)  = erfc 


Similarly,  the  covariance  of  the  motor  noise: 


(44) 


V°i(k)  = V^^(k)  + IT  p^^Var[u^(k-1 )]  i=1,...,NUM  (45) 

Both  Equations  (44)  and  (45)  require  process  (i.e.  ensemble)  statistics  at 
time  k.  However,  these  are  not  available  from  a single  Monte-Carlo 
trajectory,  and  their  precomputation  for  subsequent  read-in  is  unfeasible. 
The  approach  we  have  taken  is  to  obtain  temporal  approximations  using 
filtered  past  data.  An  approximation 
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a(k)  4 E[yj_2(k-N)] 

is  obtained  via  Ist-order  filtering  of  yj^^(k-N)  , 
a(k)  = e '^’’™a(k-D  + (1  - 

with  initial  condition  ot  ( N-1 ) = y^2(0).  The  approximate  variance  of  u^2(k-i) 
is  found  using  a two-step  procedure  that  estimates  (through  filtering)  the 
mean  and  mean-square,  and  then  computes  the  variance.  The  time  constant 

% = 0.5  sec. 

2.6.3  Pseudo-Random  Noise  Sequence 

The  Monte-Carlo  simulation  of  the  human  operator  equations  must  generate 
discrete  white-noise  sequences  for  observation  noise  and  motor  noise  that 
have  specified  variances  and  Vgj_/A  , respectively.  This  accomplished 

by  picking,  at  time  k. 


Vni 

V . (k)  = 


?(k) 


ui  - ' ^ 

where  5(k)=N(0,1)  is  a unit  variance,  zero  mean,  Gaussian  random  variable, 
(k)  is  generated  by  averaging  12  uniformly  distributed  (-1/2,  +1/2) 
independent  random  variables.  A slight  modification  is  made  when  choosing 
the  observation  noise  sequence,  to  reflect  more  precisely  the  underlying 
multiplicative  noise  process.  We  pick 

1/2 


Vyi(k)  = 


V 


yi 


^(k) 


where 


^yi(k)  = 


Vyi  +TT  Piyj_2(k-n) 
fi(k)  Ni2(k-N) 


depends  only  on  the  instantaneous  value  of  y^(k-N). 


2.7  Summary  of  Human  Model  Computations 


The  major  part  of  the  man-machine  simulation  is  the  implementation  of 
the  equations  of  the  OCM,  where  a control  input  u^^  is  generated  from  the 
observations  y^.  The  steps  in  this  process  are  summarized  below. 


12 
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1.  Storage  of  y(k)  in  the  (N+1)-st  column  of  PASTY,  Eq.  (42b). 

2.  Computation  of  a(k)  via  Eq.  (46),  and  the  observation  noise 
covariance  Vy(^)  using  Eqs . (43)-(44). 

3.  Computation  of  the  residuals 

v(k)  = y(k-N)  + ” Cp(k'k-I)  (47 

4.  Compute  the  filter  gain  via  Eq.  (3a)  and  update  the  Riccati 
equation  (40a)  to  obtain 

5.  Obtain  the  a posteriori  estimate  p(klk),  Eq.  (36a). 

6.  Obtain  the  state  estimate  x(k)  via  Eq.  (41). 

7.  Compute  the  new  commanded  control, 

u^(k)  = -L^x(k) 

and  store  it  in  the  (N+1)-st  column  of  PASTUC. 

8.  Generate  the  piecewise  constant  control  u(k)  to  use  over  the 
upcoming  time  interval  [k,k+1 ] from  Eq.  (26a)  . 

9.  Update  the  estimate  of  Var[u(k)]  for  use  at  the  next  time  step. 

10.  Propagate  E and  p using  Eqs.  (40b)  and  (36b). 

11.  Do  a stack  pushdown  (i.e.,  column  shift  to  the  left)  on  PASTY  and 
PASTUC  to  get  ready  for  the  next  time  step. 
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3.  PROGRAM  DESCRIPTION 

The.  computer  program  MCARLO  has  been  developed  for  simulating  the  human 
operator  equations,  and  controlling  the  parameter  updating  processes.  The 
program  consists  of  six  major  routines  that  are  highly  modular  in  structure: 


1. 

MAIN 

2. 

SYSTM 

3. 

UPDATE 

4. 

MAN 

5. 

INFORM 

6. 

PRINTR 

along  with  numerous  minor,  user  supplied  routines.  The  function  of  each  of 
the  above  routines  is  discussed. 

3 • 1 MAIN  Program 

The  MAIN  program  initializes  time  and  controls  the  overall  program  flow 
according  to  Figure  1.  It  calls  the  required  subroutines  for  system 
propagation,  and  information  storage.  Time  is  incremented,  t = t+A , and  the 
cycle  repeats.  When  t _>  t^^  the  printout  routine  is  called. 

3- 2 Subroutine  SYSTM 

Subroutine  SYSTM  is  a user-oriented  routine  that  simulates  the  response 
of  the  actual  system.  Given  the  control  input  u,  and  the  external  or 
disturbance  inputs  w and  z over  the' time  interval  (t-A,  t]  , SYSTM  returns  the 
value  of  y at  time  t.  At  the  first  time  step,  t = t^ ^ internal 

initializations  are  performed.  The  present  implementation  requires  SYSTM  to 
compute  and  return  (at  time  t)  the  values  of  w and  z for  its  own  use  over  the 
nejot  time  interval.  As  an  alternate  approach,  a separate  subroutine  EXTINP 
might  perform  this  function  once  NW  and  NZ  are  known. 

The  SYSTM  subroutine  is  entirely  self-contained.  As  such  it  is  possible 
to  replace  it,  in  its  entirety,  by  user  written  routine  that  simulates  the 
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given  system. (2)  The  only  requirement  is  the  generation  of  y(t)  from  control 
input  u(t)  and  the  disturbance  inputs.  If  the  system  is  time-varying,  the 
logic  for  updating  system  parameters  must  be  included  within  the  SYSTM 

routine.  The  system  simulation  can  thus  be  as  complicated  or  as  simple  as 
the  problem  may  warrant. 

The  SYSTM  subroutine  now  contained  in  the  MCARLO  program  treats  the 
general  linear  case 

x(t)  = AgX(t)  + Bgu(t)  + Esw(t)  + Fsz(t) 

x(f^)  = x(t”)  + 6x(t)  ; x(to”)  = 0 

y(t)  = C^x(t)  + Dgu(t) 

cov[W^(t)]  = W°j^(t) 

where  any  parameter  matrix  can  be  time-varying.  The  method  by  which 
parameters  are  updated  is  similar  to  the  alphanumeric  code/index  scheme  used 
in  TIVAR.  Parameters  can  be  changed  at  time  t via  external  or  card  inputs. 
They  can  also  be  changed  periodically  via  an  internal  — user  supplied  — 
subroutine  SYSNEW.  Table  la  defines  the  parameter  codes  for  SYSTM.  Note 
that  mnemonics  A and  not  AS,  etc.,  have  been  used  for  compatibility  with 
TIVAR  deck  setups. 

At  time  t the  subroutine  updates  the  pertinent  variables  and  (if 
necessary)  computes  new  discretize  equations.  The  state  is  propagated  using 
the  transition  matrix  method.  At  time  t^  matrices  A,  B,  C must  be  input  for 
proper  initialization.  Unless  input,  D,  E,  F,  Wo^  ^ assumed  to  be  zero. 
Finally,  the  parameters  associated  with  codes  1-7  are  stored  in  common 
blocks.  This  makes  them  accessible  to  other  subroutines  for  the  special  case 
when  model  = system. 


An  analog  simulation  or  a "real"  system  with  A/D  interface  provides  an 
interesting  possibility. 
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A 1 

B 2 

C 3 

D A 

E 5 

F 6 

WO  7 

XI NC  8 

INT  9 

PRINT  10 
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Table  la:  PARAMETER  CODES  IN  SYSTM 

DESCRIPTION 

System  A^  matrix,  NX  by  NX 

Control  matrix,  NX  by  NU 

Output  C^  matrix,  NY  by  NX 

Output  D matrix,  NY  by  NU 

Noise  E^  matrix,  NX  by  NW 

Bias  input  matrix  NX  by  NZ 

Noise  covariances  Wg . ^ vector 

Increment  6x  to  system  state,  NX  vector 

Transfer  to  subroutine  SYSNEW 

Printout  Interval  for  data 


CODE 

AM 

BM 

CM 

DM 

EM 

FM 

WOM 

XHINC 

TD 

MNA 

MNR 

SNA 

SNR 

TH 

ATTN 

CGAIN 

DGAIN 

INT 

PRINT 


Table  Ib:  PARAMETER  CODES  IN  UPDATE 

key  . DESCRIPTION 


1 

2 

3 

A 

5 

6 

7 

8 
9 

10 

11 

12 

13 

1 A 

15 

16 

17 

18 
19 


Model  A matrix,  NXM  by  NXM 

Control  B^  matrix,  NXM  by  NUM 

Output  C matrix,  NYM  by  NXM 

Output  raatrix,  NYM  by  NUM 

Noise  E matrix,  NXM  by  NWM 

Bias  input  raatrix,  NXM  by  NZM 

Model  noise  covariances  Wo^  vector 

Increraent  6x  to  , NXM  vector 

Huraan's  time  delay  x 

Additive  NUM  motor  noise  variances  V 

u 

Motor  noise  ratios  p^,  NUM  vector 
Additive  NYM  sensor  noise  variances  V 

y 

Sensor  ooise  ratios  ^ ^YM  vector 

Observational  thresholds,  a^^,  NYM  vector 
Attention  allocations,  f . ^ jjyM  vector 

Continuous  control  gains  , nUM  by  NTOT 
Discrete  conrol  gains  ^UM  by  NTOT 

Transfer  to  subroutine  MANNEW 
Printout  interval  for  data 
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3. 3 Subroutine  UPDATE 

The  two  major  functions  of  this  subroutine  are  to  update  the 
(time-varying)  parameters  in  the  human  operator  model,  and  to  compute  the 
discretized  model  equations.  The  parameters  are  updated  using  an 
alphanumeric  /4ode/ index  scheme . Parameters  can  be  changed  at  time  t via 
external  or  card  inputs.  They  can  also  be  changed  periodically  via  an 
internal  — user  supplied  — subroutine  MANNEW.  Table  Ib  defines  the  codes 
for  UPDATE.  The  parameters  are  described  in  Sections  2.2  - 2.3- 

The  discretization  of  the  human  operator  model  equations  follows  the 
approach  in  Sections  2.4  - 2.5-  A change  in  either  A^  gr  Bm  necessitates  a 
recoraputation  of  the  discrete  system  matrices  and/or  F^.  If  continuous 
time  feedback  gains  input,  UPDATE  computes  the  equivalent  discretized 

gains  using  the  "average  gain"  method. 

UPDATE  initializes  all  of  the  man-model  parameters  to  zero,  (3)  with  the 
exception  of  A^^^  3^^^,  Cj^,  or  which  must  be  input  at  time  t^.  The 
subroutine  can  equate  any  of  the  first  seven  code  parameters  to  their  linear 
system  counterparts. 

3.4  Subroutine  MAN 

This  is  the  major  computational  subroutine  in  the  MCARLO  program.  It 
performs  all  of  the  human  model  computations  summarized  in  Section  2.?. 

Thus,  the  basic  function  of  MAN  is  to  output  (at  time  t)  the  NU  control  u to 
apply  to  the  system  over  (t,  t+A] . The  dynamic  inputs  to  MAN  include  the  NY 
observations  y(t),  and  the  value  of  the  deterministic  input  z over  (t,  t+A]. 
This  latter  requirement  is  expected  to  be  relaxed  through  future  modeling 
efforts . 

3.5  Subroutine  INFORM 

This  subroutine  is  used  to  store  data  for  subsequent  printing  and/or 
plotting.  Data  is  stored  on  disk  files  every  (NPRNT)A seconds  starting  at 
t=t^.  For  convenience,  printed  and  plotted  variables  are  stored  on  separate 

Initial  f^  = 1 
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files.  The  value  of  NPRNT  is  an  input  parameter  to  MAIN,  but  can  be  changed 
via  either  system  or  update.  For  the  system,  any  component  of  x,  y or  u may 
be  output  as  data.  For  the  man-model,  any  component  of  x,  y or  v fnay  be 
output  where  y = C^. 

3.6  Subroutine  PRINTR 

This  routine,  called  the  first  time  t^  t outputs  the  stored  time 

r 

histories  of  the  selected  variables.  For  plotted  variables,  an  automatic 
scaling  feature  is  used. 

3.7  Program  Operation 

MCARLO  has  been  designed  to  generate  Monte-Carlo  time  histories  of  the 
signals  in  a man-machine  conrol  task.  Each  run  of  MCARLO  generates  one 
sample  path.  To  obtain  more  elements  in  the  ensemble  it  is  necessary  to  make 
additional  computer  runs,  using  different  values  for  the  random  number 
generator  seed.  The  sample  waveforms  can  all  be  stored  for  later  ensemble 
averaging.  As  the  number  of  samples  > the  sample  statistics  should 

converge  to  the  ensemble  (covariance)  statistics  that  are  computed  by  TIVAR. 
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4.  INPUT  DECK  SETUP 

There  are  three  sections  of  input  data  for  MCARLO  as  discussed  below. 

In  addition,  there  are  three  user-written  subroutines:  SYSNEW,  MANNEW,  and 

FDET. 

4. 1 Control  Cards 

There  are  5 major  control  cards  that  are  required  by  the  MAIN  program. 

Card  1 - Title  Information 
Column  1 : blank 

Columns  2-80:  alphanumeric  title  information 

Card  2 - Random  Number  Seed,  110  Format 
Field  1 : IXYZ  = any  integer 

Card  3 - Time  Information,  3^10. 0,  110  Format 

Field  1:  DEL  = discrete  time  step  (sec) 

Field  2:  TO  = initial  time  (sec) 

Field  3-  TF  = final  time  (sec) 

Field  4:  NPRNT  = printout  frequency  (integer) 

Card  4 - Print/Plot  Information  for  System  Variables,  3(2011)  Format 

Field  1 : Print/Plot  codes  for  states  1 - NX 

Field  2:  Print/Plot  codes  for  outputs  1 - NY 

Field  3:  Print/Plot  codes  for  controls  1 - NU 

Card  5 - Print/Plot  Information  for  Model  Variables,  3(2011)  Format 
Field  1 : Print/Plot  codes  for  state  estimates  1 - NXM 

Field  2:  Print/Plot  codes  for  output  estimates  1 - NYM 
Field  3-  Print/Plot  codes  for  KF  residuals  1 - NYM 

Note  that  on  cards  4-5  each  column  of  an  associated  field  corresponds  to 
one  state,  output  estimate,  etc.  A single  integer  governs  the  printing  or 
plotting  of  the  time  history  of  the  variable: 
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0,  or  blank  = no  printing  or  plotting  of  the  variable 

1 = print  time  history  vs.  time 

3 = plot  time  history  vs.  time 

2 = print  and  plot  time  history  vs.  time 

A maximum  of  10  of  any  variables  (e.g.  states  or  outputs)  can  be  printed  on 

wide  paper. 

4 . 2 System  Parameter  Cards 

These  cards  are  used  to  change  system  parameters  in  the  linear  case.  If 

the  user  supplies  his  own  SYSTM  subroutine,  this  section  of  data  is  omitted, 

or  replaced  by  problem  specific  data  cards. 

Card  1 - Frequencies  for  internal  time  breaks,  NDTS,  1015  Format 

The  10  fields  are  associated  with  the  1 Osystem  parameter  cards 
(see  Table  Ila)  on  a one-to-one  basis.  The  I-th  field  is 
associated  with  Code  I.  NDTS(I)  is  the  frequency  (number 
of  time  steps)  at  v^ich  subroutine  SYSNEW  is  to  be  called 
internally  with  KEYrI,  starting  at  time  TO.  Calling  SYSNEW 
with  KEY=I  sets  ISFLAG(I)=1  for  one  time  step.  The  actual 
parameter  values  must  be  changed  internally  by  user-written 
code.  If  no  code  is  supplied,  the  associated  parameters 
retain  their  value . 

Remaining  Cards  - These  are  used  to  change  system  parameters  via  external 

read-in  at  specified  times.  The  deck  setup  follows  a standard  form. 

Time  Card  - Cols.  1-4  Alphanumeric  TIME 

Cols.  11-20  Time  of  external  break  E10.0 

Code  Card  - Cols.  1-5  One  of  the  Alphanumeric  codes  in  Table  Ila 
Cols.  8-10  Index  NQQ  for  dimension  information,  13 

Parameter  Cards  - The  new  parameter  values  required  by  the  code. 
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Table  Ila:  SYSTEM  PARAMETER  CARD  INPUTS 


CODE 

KEY 

INDEX 

INPUT  DATA 

INITIAL  VALUE 

A 

1 

NX 

. ,NX,  j=1 , . . . ,NX 

undef 

B 

2 

NU 

(B^) . . i=1, . . 

s IJ 

.,NX,  j=1,...,NU 

undef 

C 

3 

NY 

^ ij  i= 1 • 

. ,NY,  j=1, . . . ,NX 

undef 

D 

4 

NY 

(C^) . . i=1, . . 

.,NY,  j=1,...,NU 

0 

E 

5 

NW 

^ ^s^ ij  i=  1 • 

.,NX,  j=1,...,NW 

0 

F 

6 

NZ 

(F J.  . i=1,  .. 

S IJ 

.,NX,  j=1,...,NZ 

0 

WO 

7 

NW 

(Wg)i  i=1,... 

,NW 

0 

XI NC 

8 

— 

( & ) . i=i , . . . 
1 ’ 

,NW 

0 

INT 

9 

KEY 

— 

PRINT 

10 

NPRNT  

Table  Ilb:  MAN- 

MODEL  PARAMETER  CARD  INPUTS 

CODE 

KEY 

INDEX 

INPUT  DATA 

INITIAL  VALUE 

AM 

1 

NXM* 

(A  ). . i=l, .. 
m ij 

. , NXM , j = 1 , . . . , NXM 

undef 

BM 

2 

NUM* 

(B^). . i=1,.. 
m ij 

. ,NXM,  j=1 , . . . ,NUM 

undef 

CM 

3 

NYM* 

<Vij 

. ,NYM,  j=1 , . . . ,NXM 

undef 

DM 

4 

NYM* 

<Vij 

. ,NYM,  j=1 , . . . ,NUM 

0 

EM 

5 

NWM* 

<Vij 

. ,NXM,  j=1 , . . . ,NWM 

0 

FM 

6 

NZM* 

<V1J 

. ,NZM,  j=1 , . . . ,NZM 

0 

WOM 

7 

NWM* 

(«S>i 

,NWM 

0 

XHINC 

8 

— 

i=1 , . . . 

,NXM 

0 

TD 

9 

— 

T 

0 

MNA 

10 

— 

V . i=1 , . . . , 

Ul 

NUM 

0 

MNR 

11 

— 

p . i=1 , . . . , 

Ul 

NUM  in  dB 

- dB 

SNA 

12 

— 

V . i=1,..., 
yi 

NYM 

0 

SNR 

13 

— 

Pyi  i= 1 , . . . , 

NYM 

- dB 

TH 

14 

— 

a.  i=1 , . . . ,NYM 

0 

ATTN 

15 

i=1,.. 

. ,NYM 

1 

CGAIN 

16 

(L  ). . i=1,.. 

C IJ 

. ,NUM,  j=1 , . . . ,NXM+NUM 

undef 

DGAIN 

17 

— 

. ,NUM<  j=1 , . . . ,NXM+NUM 

undef 

INT 

18 

KEY 

— 

PRINT 

19 

NPRNT 





*If  <0,  model  parameters  are  automatically  equated  to  system  parameters, 
and  no  input  data  is  needed. 
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The  sequence  of  code  card  followed  by  new  parameter  values  is  repeated 
for  all  items  that  the  user  wishes  to  change  at  the  given  time.  To  change 
parameters  at  the  next  time,  input  a new  time  card,  followed  by  a code  card, 
input  the  parameter  values,  code  card,  etc.  When  using  external  (card) 
updates,  the  following  rules  must  be  observed: 

1 . NX  + NU  •<  1 5 

2.  Time  breaks  must  occur  in  increasing  order. 

3-  Parameter  cards  should  occur  in  the  sequence  listed  in  Table  Ila . 
Thus,  codes  with  lower  KEY  numbers  should  be  read  in  first. 

4.  The  parameter  cards  must  be  input  immediately  following  the 
associated  code  card. 

5.  The  last  system  parameter  card  must  be  an  end-system  card, 
containing  the  alphanumeric  ENDS  in  cols.  1-4. 

6.  If  they  occur  at  the  same  time,  external  updates  take  precedence 
over  internal  updates. 

The  program  reads  all  of  the  system  update  cards  at  the  first  time  step, 
and  stores  all  the  information  on  a disk  file  for  sequential  read-in. 

3 Man-Model  Parameter  Cards 

These  cards  are  used  to  change  man-model  parameters  either  internally  in 
periodic  mode,  or  via  external  card  inputs.  The  deck  setup  is  virtually 
identical  in  form  to  the  previous  section . 

^ards  1-2  - Frequencies  for  internal  time  breaks,  NDTM,  1915  Format 

The  19  fields,  spread  over  2 cards,  are  associated  with  the  19 
man-model  parameter  codes  (see  Table  Ilb)  on  a one-to-one 
basis.  NDTM(I)  is  the  frequency  (number  of  time  steps)  at  which 
the  user  supplied  subroutine  MANNEW  is  to  be  called  internally 
with  KEY=I.  The  operation  is  similar  to  that  for  SYSNEW. 
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Remaining  Cards  - These  are  used  to  change  man-model  parameters  via 

external  read-in  at  specified  times.  The  deck  setup  follows 
the  standard  form  as  in  the  previous  section  4.2: 

Time  Card  - Cols.  1-4  Alphanumeric  TIME 

Cols.  11-20  Time  of  external  break  E10.0 

Code  Card  - Cols.  1-5  One  of  the  Alphanumeric  codes  in  Table  Ilb 
Cols.  8-10  Index  NQQ  for  dimension  information 

Parameter  Cards  - As  required  by  the  associated  code. 

The  sequence  of  code  card  - parameter  card  is  repeated  for  all  items  the  user 
wishes  to  change  at  the  given  time.  To  change  parameters  at  the  next  time, 
the  above  three  part  sequence  is  repeated . The  following  rules  must  be 
observed : 

1.  NXM  + NUM  15 
NUM  = NU 

NYM  = NY 

NZM  = NZ  (or  NZM  = 0) 

2.  Time  breaks  must  occur  in  increasing  order. 

3.  Parameter  codes  should  occur  in  the  sequence  listed  in  Table  lib. 
Thus,  codes  should  be  sorted  for  reading  according  to  increasing 
KEY  numbers. 

4.  The  parameter  cards  must  be  input  immediately  following  the 
associated  code  card,  unless  NQQ  < 0 for  codes  1-?. 

5.  The  last  man-model  parameter  card  must  be  an  end-man  card, 
containing  the  alphanumeric  ENDM  in  cols.  1-4. 

6.  If  they  occur  at  the  same  time,  external  updates  take  precedence 
over  internal  updates. 
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A useful  option  is  included  for  any  of  codes  1-7.  If  the  specified 
NQQ  _<  0,  man-model  parameters  are  automatically  set  equal  to  the 
corresponding  system  parameters  in  the  linear  case.  Also,  no  read-in  of 
model  parameters  is  done. 

4 . 4 Entering  Parameter  Data 

Data  is  entered  on  the  parameter  cards  in  8E10.0  Format,  i.e.,  in 
floating  point  fields  of  10  columns  with  a maximum  of  8 fields  per  card.  The 
numbers  may  be  either  in  fixed-point  (decimal)  form  or  in  scientific 
(exponential)  form  with  the  exponent  right  justified  in  the  field.  Matrices 
are  entered  one  row  at  a time.  If  a row  contains  more  than  8 entries, 
continue  on  a second  card  for  that  row.  A new  row  always  begins  on  a new 
card.  Vectors  are  entered  in  similar  8E10.0  format:  the  first  entry  in  the 

first  field,  second  entry  in  the  second  field,  etc. 

4 . 5 User  Written  Routines 

The  three  user  written  routines  are  SYSNEW(KEY) , MANNEW  (KEY),  and 
FDET(K,T).  The  purpose  of  the  first  two  routines  for  changing  system  and 
man-model  parameters  has  been  discussed  earlier.  The  function  FDET(K,T)  is 

used  to  generate  the  time  history  of  the  deterministic  (bias)  inputs 
Thus,  at  time  T,  and  for  inputs  K,  K=1,...,NZ, 

FDET(K,T)  = 

The  user  must  supply  his  own  code  for  FDET. 
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5.  SAMPLE  PROBLEM 

A sample  problem  illustrating  many  of  the  features  of  the  MCARLO  program 
is  given  in  this  section.  A description  of  the  problem  is  presented  first, 
followed  by  a listing  of  the  user  written  subroutines,  the  input  data  deck, 
and  a listing  of  the  output. 

5. 1 Sample  Problem  Description 

This  problem  analyzes  an  AAA  tracking  task.  The  controller  tracks  the 
azimuth  angle  of  a target  which  is  executing  a level  fly-by.  The  key  element 
illustrated  by  this  problem  is  the  use  of  the  FDET  function  to  generate  the 
time  history  of  the  deterministic  input,  i.e.  the  azimuth  trajectory  of  the 
target. 

The  controller  is  explicitly  presented  with  a display  of  the  azimuth 
sighting  error,  and  is  assumed  to  derive  the  corresponding  error  rate.  His 

task  is  to  minimize  this  error  by  controlling  a set  of  rate-aided 
second-order  sight  dynamics,  his  control  being  a hand-crank. 

The  target  being  tracked  is  executing  a constant  speed  straight  and 
level  fly-by  of  44  sec  duration.  The  range  of  the  target  at  crossover,  R^, 
is  3000  ft,  and  the  speed,  V,  is  733  ft/ sec  producing  a maximum  azimuth 
angular  velocity  of  about  14  deg/sec  at  crossover.  Initially,  22  sec  before 
crossover,  the  target  is  16,126  ft  from  the  crossover  point,  its  azimuth 
angular  position  is  -79.46  deg,  and  its  azimuth  angular  velocity  is  0.4683 
deg/ sec . 

The  system  states  are  defined  as  follows: 

x^  = target  azimuth  angular  position  (degrees) 

= target  azimuth  angular  velocity  (deg/ sec) 
x^  = sight  azimuth  angular  position  (degrees) 
xij  = sight  azimuth  angular  velocity  (deg/ sec) 
x^  = integral  of  the  control  input 
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It  is  assumed  that  the  controller  employs  a "constant  velocity"  model  of 

the  target  position.  Consequently,  the  state  space  equations  for  the  first 
two  states  (the  deterministic  states)  are: 


^.|(t)  = X2(t) 


x^Ct)  = z(t) , 

where  z(t),  the  deterministic  input  is  the  azimuth  angular  acceleration  of 
the  target.  Thus,  z(t)  is  given  by: 


z(t)  = -2(V/R  ',2 

C f 


D(t)/R 


where 


(1  + (D(t)/R^)2)2 


(47) 


V = the  speed  of  the  target  = 733  ft/sec 

= the  range  of  the  target  at  crossover  = 3000  ft 
D(t)  = the  distacne  of  the  target  from  the  crossover  point 

= + Vt  = -16,126  + 733t 

The  FDET  function  computes  z(t)  according  to  Eq.  47. 


The  transfer  function  relating  the  sight  position  to  the  control  input 

is: 


^2 ( _ 64( s+1 ) 

u(s)  s(s^+i 2s+64) 


(1+1/sX 


64 

s‘^+12s+64 


(48) 


Consequently,  the  state  space  equations  for  the  last  three  states  (the 
controllable  states)  are: 


X _ Y 

3 - ^4 

= -64x^(t)  - 12X|^(t)  + 64x^(t)  + 64u(t)  (49) 

*5  = u(t) 

The  displayed  outputs  are  the  azimuth  sighting  error  and  error  rate  and 
are  given  by: 
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y^(t)  = x^(t)  - x^(t) 

y^{t)  = x^ct)  - x^(t) 


(50) 


Regarding  the  human's  inherent  limitatons,  the  observation  noise  to 
signal  ratio  (SNR)  and  motor  noise  ratio  (MNR)  are  set  to  the  nominal  values 
of  -20dB  and  -25dB,  respectively.  The  perceptual  time  delay  (TD)  is  set  to 
0.20  sec.  Observational  thresholds  are  set  at  0.05  deg  for  y^  (corresponding 
to  of  the  field  of  view  of  the  gunsight) , and  0.025  deg  for 
(corresponding  to  a nominal  differential  threshold  for  motion). 


The  control  gains,  CGAIN,  computed  by  another  program,  were  chosen  by 
setting  the  cost  on  error  to  unity,  and  adjusting  the  cost  on  control  rate  to 
produce  a neuro-motor  lag  T^^  q -j  gg^  _ 

Finally,  since  the  mean  initial  states  are  far  from  zero,  the  human's 
estimator  requires  some  time  to  settle  down.  The  presence  of  the  human's 
time  delay  further  compounds  this  problem.  To  provide  an  initializing 
transition  period  while  the  human's  estimator  settles  down,  the  sample 
problem  is  started  at  t=-6.0  sec,  at  which  time  the  target  is  20,524  ft  from 
the  crossover  point,  its  azimuth  angular  position  is  -81.684  deg  and  its 

azimuth  angular  velocity  is  0.2928  deg/ sec.  During  this  initialization 
period,  the  human's  time  delay  is  set  to  zero,  and  printouts  are  suppressed. 
At  t=0,  however,  the  time  delay  is  set  to  0.20  sec,  and  the  printout  interval 
is  set  to  1 sec. 

The  states  are  incremented,  XING,  so  that  the  initial  angular  position 
and  velocity  of  the  target  are  correct,  and  the  human's  state  estimates  are 
incremented,  XHINC,  so  that  the  initial  error  and  error  rate  are  zero. 
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5 . 2 User  Written  Subroutines  for  the  Sample  Problem 

MCP  - PROBLEM  DEPENDENT  SUBPROGRAMS  FOR  MCARLO 
INCLUDES: 

1 - FUNCTION  FDET 

2 - SUBROUTINE  SYSNEW 

3 - SUBROUTINE  MANNEW 


FUNCTION  FDET(NQT,T) 

C GENERATES  DETERMINISTIC  INPUT 

DATA 

1 XO,  YO,  VO  /-16126.0,  3000.0,  733-0/, 

2 R 757.296/ 

X=X0+V0*T 

A=X/Y0 

B=1 .0+A»A 

C=(V0/Y0)/B 

D=-2.0»A»C»C 

FDET=R»D 

RETURN 

END 


SUBROUTINE  SYSNEW(NQQ) 

C INTERNAL  UPDATES  TO  THE  SYSTEM 

RETURN 

END 


SUBROUTINE  MANNEW(NQQ) 

C INTERNAL  UPDATES  TO  THE  MAN 

RETURN 

END 


- 28 


Report  No.  3463 


Bolt  Beranek  and  Newman  Inc. 


5.3  Input  deck  for  the  Sample  Problem 

P-I-D  CONTROLLER.  MONTE-CARLO  SIMULATION 
56315 
0.05 


22221 

10100 

0 


0 

5 

0.0 

0.0 

0.0 

0.0 

0.0 

1 

0.0 

0.0 

0.0 

64.0 

1.0 

1 

0.0 
1 .0 
0.0 
0.0 
0.0 
2 

1.0 

0.0 

2 

0.0 

0.0 


C 
D 

XI NC 

-81.684 
ENDS 

0 
0 

AM 

BM 
CM 
DM 
FM 

™§1.684 

TD 

0.00 

CGAIN 

-15.39 

10.00 


0 

0 

0 


0 

0 


-6, 

,0 

44, 

,0 

120 

22 

2 

20 

20 

0 

0 

0 

0 

0 

0 

0 0 

1 . 

.0 

0. 

.0 

0.0 

0, 

,0 

0. 

.0 

0. 

.0 

0.0 

0. 

,0 

0 

.0 

0, 

.0 

1.0 

0. 

.0 

0 

.0 

-64, 

,0 

-12.0 

64, 

.0 

0 

.0 

0 

,0 

0.0 

0, 

.0 

0.0 

1.0 


0.2928 


0 

0 


0.2928 


■3.303 


MNR 

•25.0 

SNR 

•20.0 

-20.0 

TH 

0.05 

0.025 

TIME 

0 

0.0 

TD 

0.20 

PRINT 

ENDM 

20 

-1 .0 

0.0 


-81.684 
0 0 


0.0 

-1.0 


0.2928 
0 0 


-81.684  0.2928 


6.359  0.6404 


0.0 

0.0 


0.0 

0 


0.0 


9.034 
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5-4  Output  listing  for  the  Sample  Problem 

STARTING  MCARLO 
2- Dec-76  14:20 


P-I-D  CONTROLLER.  MONTE-CARLO  SIMULATION 


RANDOM  NUMBER  SEED=  56315 

INTEGRATION  TIME  STEP  = 0.050 

INITIAL  TIME  = -6.000 

TERMINAL  TIME  = 44.000 

PRINTOUT  FREQUENCY  = 120 

SYSTEM  INTERNAL  BREAKS  INDEX  CODE 


NDT 


SYSTEM  EXTERNAL 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 


BREAK  AT  T= 
1 .OOOE+00 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 


-6.000 

O.OOOE-01 

O.OOOE-01 

O.OOOE-01 

-6.400E+01 

O.OOOE-01 


CODE  A 

O.OOOE-01 
O.OOOE-01 
1 .OOOE+00 
-1 .200E+01 
O.OOOE-01 


INDEX=  5 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 
6.400E+01 
O.OOOE-01 


SYSTEM  EXTERNAL  BREAK 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 
6.400E+01 
1 .OOOE+00 


AT  T=  -6.000  CODE 


B INDEX=  1 


SYSTEM  EXTERNAL  BREAK 
O.OOOE-01 
1 .OOOE+00 
O.OOOE-01 
O.OOOE-01 
O.OOOE-01 


AT  T=  -6.000  CODE 


F INDEX=  1 


SYSTEM  EXTERNAL  BREAK  AT  T= 
1.000E+00  O.OOOE-01 

O.OOOE-01  1.000E+00 


-6.000  CODE  C 
-1.000E+00  O.OOOE-01 
O.OOOE-01.  -1.000E+00 


INDEX=  2 
O.OOOE-01 
O.OOOE-01 


SYSTEM  EXTERNAL  BREAK  AT  T=  -6.000  CODE  D INDEX-  2 

O.OOOE-01 
O.OOOE-01 


SYSTEM  EXTERNAL  BREAK  AT  T= 
-8.168E+01  2.928E-01 

HUMAN  INTERNAL  BREAKS  INDEX 


-6.00C 

-8.16J 

CODE 


CODE 

E+01 

NDT 


2. 


XI NC 
928E-01 


INDEX=  0 
O.OOOE-01 


HUMAN  EXTERNAL  BREAK 
HUMAN  MODEL  = SYSTEM 

HUMAN  EXTERNAL  BREAK 
HUMAN  MODEL  = SYSTEM 

HUMAN  EXTERNAL  BREAK 
HUMAN  MODEL  = SYSTEM 

HUMAN  EXTERNAL  BREAK 
HUMAN  MODEL  = SYSTEM 

HUMAN  EXTERNAL  BREAK 
HUMAN  MODEL  = SYSTEM 

HUMAN  EXTERNAL  BREAK 
-8.168E+01  2.9; 


AT 

T= 

-6.000 

CODE 

AT 

T= 

-6.000 

CODE 

AT 

T= 

-6.000 

CODE 

AT 

T= 

-6.000 

CODE 

AT 

T= 

-6.000 

CODE 

AT 

8E- 

T= 

-01 

-6.000  CODE 

-8. 168E+01 

AM 

INDEX= 

0 

BM 

INDEX= 

0 

CM 

INDEX= 

0 

DM 

INDEX= 

0 

FM 

INDEXr 

0 

XHINC 

INDEX= 

0 

2.928E-OI  O.OOOE-01 
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HUMAN  EXTERNAL 
O.OOOE-01 

BREAK  AT  T= 

-6.000 

CODE 

TD 

INDEX= 

0 

HUMAN  EXTERNAL 
-1 .539E+01 
1 .OOOE+01 

BREAK  AT  T= 
-3.3O3E+OO 

-6.000  CODE 

6.359E+00 

CGAIN 

6.404E-01 

INDEX=  0 

9.034E+00 

HUMAN  EXTERNAL 
-2.5OOE+OI 

BREAK  AT  T= 

-6.000 

CODE 

MNR 

INDEX= 

0 

HUMAN  EXTERNAL 
-2.000E+01 

BREAK  AT  T= 
-2.000E+01 

-6.000 

CODE 

SNR 

INDEX= 

0 

HUMAN  EXTERNAL 

BREAK  AT  T= 

-6.000 

CODE 

TH 

INDEXz 

0 

5.000E-02  2.500E-02 

EQUIVALENT  DISCRETE  GAINS  GENERATED: 
-1.187E+01  -2.871E+00  " 

8.741E+00 


HUMAN  EXTERNAL  BREAK  AT  T= 
2.000E-01 


4.078E+00 
0.000  CODE 


4.587E-01 

TD 


7.79IE+OO 
INDEX=  0 


HUMAN  EXTERNAL  BREAK  AT  T= 
MEAN  OF  X VECTOR 

-9.974E+00  3.223E+00 

RMS  VALUES  OF  X VECTOR 

6.586E+01  3.843E+00 

MEAN  OF  Y VECTOR 

-1.112E-01  2.335E-02 

RMS  VALUES  OF  Y VECTOR 

2.074E+00  1.195E+01 

MEAN  OF  U VECTOR 
1.580E+00 

RMS  VALUES  OF  U VECTOR 
9.259E+00 


0.000  CODE  PRINT  INDEX=  20 
-9.863E+OO  3-199E+00  -I.O8OE+OI 

6.585E+01  1.257E+01  6.407E+01 
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TIME  STATE  1 

0.00  -7.947E+OI 

1-00  -7.898E+01 

2.00  -7.844E+01 

3.00  -7.786E+OI 

4.00  -7.72IE+OI 

5.00  -7.648E+01 

6.00  -7.567E+01 

7.00  -7.477E+01 

8.00  -7.3711E+OI 

9.00  -7.257E+01 

10.00  -7.122E+01 

11.00  -6.965E+OI 

12.00  -6.782E+OI 

13.00  -6.564E+01 

14.00  -6.3O3E+OI 

15.00  -5.984E+01 

16.00  -5.59OE+OI 

17.00  -5.095E+01 

18.00  -4.467E+01 

19.00  -3.667E+OI 

20.00  -2.658E+OI 

21.00  -I.436E+OI 

22.00  -6.721E^01 

23.00  I.3IOE+OI 

24.00  2.55OE+OI 

25.00  3.58IE+OI 

26.00  4.401E+01 

27.00  5.045E+01 

28.00  5.55IE+OI 

29.00  5.954E+01 

30.00  6.279E+01 

31.00  6.546E+01 

32.00  6.767E+01 

33.00  6.954E+01 

34.00  7.II3E+OI 

35.00  7.25OE+OI 

36.00  7.368E+OI 

37.00  7.472E+OI 

38.00  7.564E+01 

39.00  7.646E+01 

40.00  7.719E+01 

41.00  7.785E+01 

42.00  7.844E+01 

43.00  7.898E+OI 


STATE  2 
4.677E-OI 
5.116E-01 
5.6I8E-OI 
6. I97E-OI 

6.868e-01 

.653E-01 
.577E-01 
9.675E-01 
1 .099E+00 
1 .258E+OO 
1 .454E+00 
1 .696E+OO 
2.000E+00 
2.388E+OO 
2.89OE+OO 
3.548E+OO 
4.421E+00 
5.584E+00 
7.117E+00 
9.054E+00 
125E+01 
317E+OI 
400E+01 
325E+01 
I36E+OI 
. - I6OE+OO 
7.204E+00 
5.65IE+OO 
4.471E+00 
3.586E+OO 
2.919E+00 
2.410E+00 
2.018E+00 
1 .71 OE+00 
1 .465E+00 
1 .267E+00 
1 . 106E+00 
9.735E-OI 
8.627E-01 
7.696E-OI 
6.905E-01 
6.228E-01 
5.645E-01 
5.139E-01 


STATE  3 
-7.96OE+OI 
-7.915E+01 
-7.845E+01 
-7.789E+01 
-7.726E+01 
-7.651E+01 
-7.565E+01 
-7.480E+01 
-7.383E+01 
-7.258E+OI 
-7. 121E+01 
-6.964E+01 
-6.78OE+OI 
-6.563E+01 
-6.297E+01 
-6.001E+01 
-5.59IE+OI 
-5. 1 14E+01 
-4.506E+01 
-3.69OE+OI 
-2.679E+01 
-1 .51 1E+01 
-1 . 157E+00 
1 .343E+OI 
2.582E+01 

3.6O7E+OI 

4.423E+OI 

5.1 16E+01 
5.6O3E+OI 
5.972E+OI 
6.31 1E+01 
6.549E+01 
6.78IE+OI 
6.959E+01 
7.I2IE+OI 
7.265E+01 


374E+01 

479E+OI 

,569E+01 

.647E+01 

.721E+01 

787E+01 

841E+01 

896E+01 


STATE  4 
4.125E-01 
>71E-01  -7 

547E-01  -7 

7.59IE-OI 
7.768E-OI 
8.245E-01 
8.815E-01 
7.255E-01 
1 .253E+00 
1 .210E+00 
1 .464E+00 
1 .852E+00 
1 .842E+00 
2.699E+00 
2.853E+OO 
4.086E+00 
4.59OE+OO 
5.437E+OO 
7.672E+00 
8.665E+00 
1 . 159E+01 
1 .323E+OI 
1 .569E+01 
1 .394E+01 
1 .O36E+OI 
8.595E+00 
6.746E+00 
4.902E+00 
4.033E+00 
3.544E+00 
3.58OE+OO 
2.432E+00 
1 .735E+00 
1 .648E+00 
1 .437E+OO 
1 .688E+00 
1 . 194E+00 
8.909E-01 
8.773E-01 
7.686E-OI 
6.342E-01 
325E-01 


STATE  5 
7.974E+01 
3E+0 1 
, . --8E+01 
-7.833E+01 
-7.773E+01 

-7.707E+OI 

-7.629E+OI 

-7.546E+01 

-7.458E+OI 
-7.349E+OI 


-7 

-7 

-6 

-6 

-6 

-6 


^.'§7^1-01 


-27E+01 
086E+01 

92IE+OI 

729E+01 
499E+01 
--228E+01 
-5.887E+01 
-5.472E+01 
-4.956E+01 
-4.282E+01 
-3.428E+OI 
-2.38OE+OI 

-1 . 150E+01 
2. 197E+00 
1 .535E+01 

2.7O8E+OI 
3.684E+01 
4.485E+01 
5. IO8E+OI 
5.599E+01 
6.001E+01 
6.314E+01 
6.58OE+OI 
6.796E+OI 
6.979E+01 
7. I36E+OI 
7.27OE+OI 
7.388E+OI 
7.489E+01 
7.579E+01 
7.659E+01 
7.73IE+OI 
7.793E+01 
7.85OE+OI 
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TIME 

0.00 

1.00 

2.00 

3.00 

4.00 

5.00 

5.00 

7.00 
8.00 
9.00 

10.00 

11.00 

12.00 

13.00 

14.00 

15.00 
16.00 

17.00 

18.00 
19.00 
20.00 
21.00 
22.00 
23.00 

24.00 

25.00 

25.00 

27.00 

28.00 

29.00 

30.00 

31.00 

32.00 

33-00 

34.00 

35.00 

36.00 

39.00 

40.00 

41.00 

42.00 

43.00 


OUTPUT  1 
1 .374E-01 
1 .656E-01 

2.328E-O3 

3.O37E-O2 
5.052E-02 
2-552E-02 
-2.782E-02 
3.871E-02 
8.75OE-02 
1 . 187E-02 

-3.47OE-O3 
-1 .513E-O2 
-1 .55IE-02 

-1 . 1 12E-02 
-5.234E-02 
1 .71  IE-01 
1 .5I3E-O2 
1 .836E-01 
3.834E-01 
2.285E-01 
2.084E-01 
7.453E-01 
4.852E-01 
-3.329E-OI 
-3. 137E-01 
-2.618E-01 
-2.212E-01 
-7.1 17E-01 
-5. 173E-01 
-1 .845E-01 
-3. 179E-01 
-3.509E-O2 
-1 .338E-OI 
-5.40QE-Q2 
-7.9I7E-O2 
-1 .501E-01 
-5.9OIE-O2 
-7.014E-02 
-4. 147E-02 
-1 .335E-02 
-1 .33OE-O2 
-1 .702E-02 
3.859E-02 

2.873E-O2 


OUTPUT  2 
5.519E-02 
-4. I55E-OI 
-2.293E-02 
-1 .395E-01 
-9.OOIE-O2 
■02 
■02 
■01 

-1 .538E-01 
4.863E-02 
-9.822E-03 
-1 .561E-01 
1 .587E-OI 
-3. 106E-01 
3.662E-02 
-5.383E-OI 
-1 .693E-01 
1 .468E-01 
-5.552E-01 
3.886E-OI 
-3.483E-01 
-5. 167E-02 
-1 .686E+00 
-6.948E-01 
9.963E-OI 
5.649E-01 
4.583E-01 
7.49IE-OI 
4.38OE-OI 
4.247E-02 
-6.614E-01 
-2. I3OE-O2 
2.827E-01 
6.138E-Q2 
2.728E-O2 
-4.204E-01 
-8.793E-O2 
8.256E-02 
-1 .45bE-02 
9.882E-04 
5.629E-02 

i.mul 

2.686E-02 
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TIME  CONTROL  1 
0.00  1.485E-01 

1.00  4.880E-01 

2.00  5.266E-01 

3.00  6.O98E-OI 

4.00  6.487E-01 

5.00  7.193E-01 

6.00  6.189E-01 

7.00  8.221E-01 

8.00  I.OI6E+OO 

9.00  I.137E+OO 

10.00  1.303E+00 

11.00  1.553E+00 

12.00  1.752E+00 

13.00  2.160E+00 

14.00  2.563E+00 

15.00  3.119E+00 

16.00  3.814E+00 

17.00  4.742E+00 

16.00  5.922E+00 

19.00  7.472E+OO 

20.00  9.74OE+OO 

21.00  1.151E+01 

22.00  1.352E+01 

23.00  1.364E+01 

24.00  1.246E+01 

25.00  1.050E+01 

26.00  8.561E+00 

27.00  6.933E+00 

28.00  5.32IE+OO 

29.00  4.479E+00 

30.00  3-609E+00 

31.00  2.83IE+OO 

32.00  2.292E+00 


35.00 

36.00 

37-00 

36.00 

39.00 

40.00 

41.00 

42.00 

43.00 


1 .939E+QQ 
1 .693E+00 
1 .536E+00 
1 .241E+00 
1 .073E+OO 
9.634E-01 
6.-557E-01 
7.456E-OI 
6.433E-OI 
5.878E-OI 
5.449E-01 
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9.00  -7 

0.00  -7 


TIME 

0.00 

1.00 

2.00 

3.00 

4.00 

§:88 

7.00 
8.00 
. .00 

10.00 

11.00 

12.00 

13.00 

14.00 

15.00 

16.00 

17.00 

18.00 

19.00 

20.00 
21.00 
22.00 

23.00 

24.00 

25.00 
26.00 

27.00 

28.00 

29.00 

30.00 

31.00 

32.00 

33.00 

34.00 

35.00 

36.00 

^^:88 

39.00 

40.00 

11:88 

43.00 


XMHAT  1 
-7.975E+01 
■7.91  1E+01 
-7.847E+01 
■7.782E+01 
-7.7I8E+OI 


7.478E+OI 
7.372E+01 
.256E+OI 
. 122E+01 
6.962E+01 
6.784E+01 
6.56OE+OI 
6.290E+01 
■5.98IE+OI 
■5.584E+01 

:^:8§fE:8] 

-3.670E+01 

■2.652E+01 
■1 .464E+01 
-6.236E-01 
1 .317E+OI 
2.546E+01 
3.575E+OI 
4.408E+01 
5.O2OE+OI 
5.542E+01 
5.973E+OI 
6.3O8E+OI 
6.570E+01 
6.786E+OI 
6.973E+OI 
7. I3OE+OI 
7.276E+OI 
7.385E+OI 

7.66OE+OI 

7.73IE+OI 

?:?§p:81 

7.9O6E+OI 


:? 


XMHAT  3 
■7.98OE+OI 
-7.924E+01 
-7.848E+01 
-7.786E+OI 
■7.723E+OI 

:?:88?i:81 

■7.477E+OI 
7.38OE+OI 
.257E+01 
.I26E+OI 
6.965E+01 
6.782E+01 
6.568E+OI 
-6.296E+OI 
-5.998E+OI 
-5.59OE+OI 

;5:8?i:81 

-3.685E+OI 
-2.672E+01 
-1 .523E+01 
-1 .045E+00 
1 .346E+01 
2.585E+OI 
3.6IOE+OI 
4.419E+01 
5.1 10E+01 
5.599E+01 
5.994E+01 
6.326E+01 
6.571E+01 
6.801E+01 
6.974E+01 
7. 137E+01 
7.279E+01 
7.387E+OI 

HWM 

7.657E+01 

7.73IE+OI 

?;581i:81 

7.9O6E+OI 
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TIME 

Y 

0. 

00 

5. 

1 . 

00 

1. 

2. 

00 

6 . 

3. 

00 

3. 

4. 

00 

4. 

5 . 

QQ 

1 . 

6 . 

00 

1 . 

7. 

00 

-4. 

8. 

00 

7. 

9. 

00 

1 . 

10. 

00 

-2. 

1 1 . 

00 

2. 

12. 

00 

-2. 

1 3 • 

00 

7. 

1 4 . 

00 

6 . 

1 5 . 

00 

1 . 

1 6 . 

00 

5. 

00 

1 . 

18. 

00 

3. 

19. 

00 

1 . 

20. 

00 

2. 

21 . 

00 

5. 

22. 

00 

4. 

23. 

00 

-2. 

24. 

00 

-3. 

25. 

00 

-3. 

26. 

00 

-1  . 

21- 

00 

-9. 

28. 

00 

-5. 

29- 

00 

-2. 

30. 

00 

-1  . 

31  . 

00 

-1  . 

32. 

00 

-1  . 

33- 

00 

_3 . 

34. 

00 

-6 . 

35. 

00 

-2. 

36. 

00 

-2. 

3Z- 

00 

_3 , 

jO  • 

00 

3 . 

39. 

00 

2. 

40. 

00 

1 . 

41 . 

00 

-2. 

42. 

00 

1 . 

43. 

00 

1 . 

MHAT  1 
440E-02 
286E-01 
689E-03 
072E-02 
954E-02 
996E-02 
385E-02 
005E-03 
715E-02 
180E-02 
016E-02 
993E-02 


40 


E-02 
OOE-02 
OQE-01 
68E-02 

98E-01 

62E-01 


901E-01 

217E-01 

8fa1E-01 


552E-02 

482E-01 

SI3E-O3 

8O6E-O2 

637E-02 

666E-02 

475E-02 

157E-02 

213E-O2 

IOIE-O3 

96IE-O2 

lOOE-02 
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TIME 

0.00 

1.00 


2.00 

3.00 

4.00 

5.00 

6.00 

7.00 

8.00 

9.00 
10.00 
11.00 
12.00 
13-00 

14.00 

15.00 

16.00 


17.00 

18.00 


00 
.00 
21.00 
22.00 

23.00 

24.00 

25.00 

26.00 

27.00 

28.00 

29.00 

30.00 

31.00 

32.00 

33.00 

34.00 

35.00 

3b. 00 


37.00 

38.00 

39.00 

40.00 

41.00 

42.00 

43.00 


INOVAT  1 
4.335E+03 
8.887E-02 
-2.526E+02 
-6.851E-02 
-2.569E-01 


1 .35OE-OI 
2.387E-02 
-7.826E-OI 
2.523E-01 
1 . 157E+03 
-9.278E+00 
3.482E-01 
2.497E-02 
-2.527E-01 
8.5O3E+O2 
-1 .96OE-OI 
2.483E-01 
-3.536E-02 


3.922E-OI 
1 .56IE-OI 
-1 .411E-01 
6.618E-01 
2.696E-01 
-5.703E-O2 
-7.549E-OI 
-4. I78E-OI 


2.654E-01 
-1 .585E-01 
2.572E-01 
-1 .941E-01 


\2OE-O2 
-9.707E-02 

-1 .796E-01 
1 . 1 17E-01 
-2.200E-01 


-5.386E-O2 
-5.22QE-02 
1 .576E-01 
3.34OE-OI 
-9.23OE-O2 


E-01 

E+00 
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STATE  1 


- 38  - 


Report  No.  3^63 


Bolt  Beranek  and  Nevman  Inc. 


X O.OOE-01 


X 4.30E+01 

STATE 


4.677E-01  X X 1.400E+01 


2 
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X O.OOE-01 


X 4.30E+01 
STATE 
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X O.OOE-01 


X 4.30E+01 
STATE 


4.125E-01  X 


X 1.569E+01 


+ X 
+X 
+X 
+ 

+X 
+ X 
+ X 

+x 


+ 

+ 

+ 


X 

X 


X 

X 


+ 

+ 


+ 

+ 


+ 

+ 


+ 

+ 


+ 

+ 

+ 


+ 

+ 

+ 

+ 

+ 

+ 

+ 


+ 

+ 

+ 


+ 

+ 

+ 


+ 

+ 

+ 

+ 

+ 

+ 


+ 

+ 

+ 


+ 

+ 

+ 


+ 

+ 

+ 

+ 

+ 

X+ 

+ 

+ 


X 

+ 


+ 

+ 


X 

X 


X 

X 


+ 

+ 

+ 


X 

X 


+x 


+ 

+ 

+ 


+ 

+ 

+ 


+x 

X 

X 


+ 

+ 
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X O.OOE-01 


X 4.30E+01 

OUTPUT 


-7.117E-01  X X 7.453E-01 


1 


42 
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X O.OOE-01 


X 4.30E+01 

OUTPUT 


-1.686E+00  X X 9.963E-OI 


2 
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X O.OOE-01 


X 4.30E+01 

CONTROL 


1.485E-01  X 


X 1.364E+01 
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X 


X 


O.OOE-01 


4.3OE+OI 

YMHAT 


-9.002E-01  X 

X 5.9OIE-OI 

+ ++++++++++++++++ 

+ 

+ 

+x 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

X 

, + 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

x+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

x+ 

+ 

+ 

+ 

x+ 

+ 

+ 

+ 

+ X 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

x+ 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ X 

+ 

+ 

+ 

+ X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ X 

+ 

+ 

+ 

+x 

+ 

+ 

+ 

+ X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

X 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ X 

+ 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

> 

: + 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

+ 

+ 

X 

+ 

+ 

1 
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X O.OOE-01 


-2.526E+02  X 


+ 

X 


X 

X 

X 

X 

X 

X 

X 


X 4.335E+03 


+ 

+ 


X 

X 

X 

X 

X 

X 

X 


+ 

+ 

+ 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


X 4.3OE+OI 


X 

X 

X 

X 

X 

X 

X 

X 


INOVAT  1 
STOPPING  MCARLO 
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6.  COMMON  BLOCK  USAGE 

Named  COMMON  blocks  are  used  to  store  most  data  arrays  and  to  pass 
information  among  the  various  subroutines.  These  are  described  below. 

1.  /PL0T1/ 

Required  inputs  for  lineprint  plot  subroutine 

2.  /INOU/  KIN,  KOUT,  KPTR,  KPUNCH,  KDISK,  IPOS,  IGOS 

Logical  unit  numbers  for  I-O  devices  and  disks  for  storage  of 
system  parameters  and  generated  data. 

3.  /INFO/ 

Storage  of  print/plot  information  including  number  of  variables, 
min  and  max  values  for  plot  scaling,  etc. 

4.  /TIMES/  TIME,  DEL,  TO,  TF , NPRNT 

Time  information  t.  A,  t^,  t^,  printout  frequency 

5.  /MAIN!/  NDIM,  NDIMl , COM!  /MAIN2/  COM2 

Common  blocks  required  for  library  subroutines 

6.  /SYSX/  NX,  NU,  BS,  AS 

Linear  System  state  parameters  BS=B^,  AS=Ag 

7.  /SYS AD/  BD,  AD 

Discrete  System  parameters  BD=r^^  AD=^g 

8.  /SYSY/  NY,  DS,  CS 

Linear  System  output  parameters  DS=D  CS=C 

s ’ s 

9.  /SYSW/  NW,  WO,  ES  /SYSZ/  NZ , FS 

External  System  input  parameters  W0=Wo  eS=E  , FS=F 

s’  s s 

10.  /SYS INC/  XI NC 

State  increment  6x 

11.  /MANX/  NXM,  NUM,  BM,  AM 

Man-model  state  parameters  BM=B  , AM=A 

' m ’ m 
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12.  /MANAD/  BD,  AD 

Augmented  discretized  Man-model  parameters  BD=F , AD=f 

13.  /MANY/  NYM,  DM,  CM 

Man-model  output  parameters  DM=D  , CM=C  (augmented) 

m’  m ^ 

in.  /MANW/  NWM,  WOM,  EM  /MANZ/  NZM,  FM 

Man-model  external  input  parameters  W0M=W°,  EM=E  , FM=F 

“>  m m 

15.  /RATIOS/  PU,  VU,  PY,  VY,  TH,  ATTN,  SIGMA 

Model  parameters  p^,  , p^,  V^,  a,  f, 

16.  /MANINC/  TD,  NPRED,  XHINC 

Man-model  parameters  TD=t , NPRED=[x/A ]=N , XHINC=6p 

17.  /GAINBK/  CGN 

CGN  = discrete  control  gains  or  equivalent  discrete 
gains  computed  from 


- ns  - 
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Report  No.  3^63 


Bolt  Beranek  and  Newman  Inc. 


C NO  TABS 


7.  MCARLO  LISTING 


MCARLO  - TIME  VARYING  MONTE  CARLO  MAN/MACHINE  SIMULATION 


INCLUDES 

1 - BLOCK  DATA  MCDAT 

2 - MAIN. 

3 - SUBROUTINE  MCARLO 


- INITIALIZES  VARIOUS  COMMON  BLOCKS 

- CALLS  SUBROUTINE  MCARLO 

- PRIMARY  SUBPROGRAM 


ALSO  REQUIRES  THE  FOLLOWING  SUBPROGRAM  FILES 


MCCMP  - COMPUTATIONS  FOR  MCARLO 
INCLUDES 

1 - SUBROUTINE  SYSTM  - PROPAGATES  THE  SYSTEM  S RESPONSE 

2 - SUBROUTINE  MAN  - PROPAGATES  THE  MAN'S  RESPONSE 

3 - FUNCTION  GAUSS  - PICKS  A NUMBER  FROM  A GAUSSIAN 

DISTRIBUTION 


MCP  - PROBLEM  DEPENDENT  SUBPROGRAMS  FOR  MCARLO 
INCLUDES 

1 - SUBROUTINE  FDET  - GENERATES  DETERMINISTIC  INPUT 

2 - SUBROUTINE  SYSNEW  - INTERNAL  UPDATES  TO  THE  SYSTEM 

3 - SUBROUTINE  MANNEW  - INTERNAL  UPDATES  TO  THE  MAN 


MCIO  - I/O  FOR  MCARLO 
INCLUDES 

1 - SUBROUTINE  UPDATE 

2 - SUBROUTINE  INFORM 

3 - SUBROUTINE  PUTOUT 

4 - SUBROUTINE  PRINTR 


- PERFORM  EXTERNAL  UPDATES 

- DO  OUTPUT  FOR  A SINGLE  TIME  STEP 

- SAVES  OUTPUT  ON  FILES 

- PRINT  THE  OUTPUT  AT  THE  END  OF  A RUN 


KPLOT  - LINEPRINTER  PLOTTING  PACKAGE 
INCLUDES 

1 - SUBROUTINE  KPLOT 

2 - SUBROUTINE  ADJUST 

3 - SUBROUTINE  QINIT 

4 - SUBROUTINE  KPLOTC 

5 - SUBROUTINE  QPLOT 

6 - SUBROUTINE  PLACE 

7 - SUBROUTINE  QPRINT 

BLOCK  DATA  MCDAT 
C INITIALIZES  VARIOUS  COMMON  BLOCKS 

COMMON 

1 /INOU/  KIN,  KOUT,  KPTR,  KPUNCH, 

1 KDIoK  IPOS  IGOS 

2 /PLOT1/  NV,  NH,  NCPW,  LW,  XL,  XH,  YL,  YH,  NXES,  NDIR,  1ST, 

2 NGLV,  NGLH,  BSYM,  GSYM,  PSYM,  ND1 , ND2 , NOUT 

DATA 

1 NV,  NH,  NXES,  NDIR,  NGLV,  NGLH,  BSYM,  GSYM,  PSYM,  1ST,  ND1 
1 / 51,  101,  1,  10,  20,  20,  1H+,  1H+,  1HX,  11,  1/ 
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PROGRAM  MCMAIN 

1 (INPUT,  OUTPUT,  TAPE5=INPUT,  TAPE6=0UTPUT, 

2 DISK,  POS,  GOS,  TAPE7=DISK,  TAPE8=P0S,  TAPE9=G0S) 

MAIN  PROGRAM 

CALLS  SUBROUTINE  MCARLO 


CALL  MCARLO 
END 
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SUBROUTINE  MCARLO 
PRIMARY  SUBPROGRAM 

DIMENSION 

1 TITLE(8), 

2 X(30),  y(30),  U(30).  W(30),  Z(30), 

3 XMHAT(30),  YMHAT(30),  RES(30) 


1 

1 

2 

2 

3 

4 


6 

7 

7 

Z 


COMMON 

/INOU/  KIN,  KOUT,  KPTR,  KPUNCH, 

KDISK  IPOS  IGOS 

/PLOT1/  NV,  NH,  NCPW,  LW,  XL.  XH,  YL,  YH.  NXES,  NDIR, 
NGLV,  NGLH,  BSYM,  GSYM,  PSYM,  ND1 , ND2,  NOUT 
/MAIN1/  NDIM,  NDIM1,  COM1(15,15) 

/MAIN2/  C0M2(15,15),  ST0RE(1800) 

/TIMES/  TIME,  DEL,  TO,  TEND,  NPRNT  ^ ^ ^ 

/INFO/  NREC,  NPRINT,  NPLOT.  LPRNTS(60),  LPRNTM(60), 
IP(6),  IG(6),  SMIN(^I),  SMAX(21) 

/FILES/  KKB,  NAMIN,  NAMOUT,  NAMDYN,  NAMPCH,  NAMDSK, 
NAMPOS,  NAMGOS 

/COMP1/  X,  Y,  U,  W,  Z,  XMHAT,  YMHAT , RES 


DATA 

1 KIN,  KOUT,  KPTR,  NOUT  /5,  6,  6,  6/, 

2 KPUNCH,  KDISK,  IPOS,  IGOS  /6 , 7,  8,  9/ 


C SET  NDIM 

NDIM=15 
NDIM1=NDIM+1 

C WRITE  DAYTIM 

150  CALL  PAGEFD(KOUT.I) 

WRITE  (KOUT, 1500) 

1500  FORMAT  (1H  ,/,1H  ,15HSTARTING  MCARLO) 

CALL  DAYTIM(KOUT) 

C ZERO  THE  VECTORS 

200  DO  220  1=1,30 

X(I)=0.0 
Y(I)=0.0 
W(I)=0.0 
Z(I)=0.0 
U(I)=0.0 
220  CONTINUE 

ISTEP=0 


C 

300 

1040 

350 

1045 

1051 

1052 


GET  TITLE  AND  RANDOM  NUMBER  SEED 
READ  (KIN, 1040)  (TITLE(l),  1=1,8) 

IF  (EOF(KIN))  1000,350 
FORMAT  (8A10) 

WRITE  (KOUT, 1045)  (TITLE(I),  1=1,8) 
FORMAT  (/,1H  ,8A10,/,1H  ) 

READ  (KIN,  1051)  IDUM 
FORMAT  (110) 

WRITE  (KOUT, 1052)  IDUM 

FORMAT  (21H  RANDOM  NUMBER  SEED=  ,110) 

CALL  RANSET(IDUM) 


C SPECIFY  DEL,  TO,  TEND.  AND  NPRNT 

330  READ  (KIN. 1060)  DEL,  TO,  TEND,  NPRNT 

1060  FORMAT  (3E10. 0,110) 

TEND=IFIX( ( TEND-TO+0 . 000 1 ) /DEL) *DEL+TO 
WRITE  (KOUT, 1065)  DEL.  TO,  TEND,  NPRNT 
1065  FORMAT  (25H  INTEGRATION  TIME  STEP  = ,F10.3,/, 
1 17H  INITIAL  TIME  = ,F10.3,/, 
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2 
3 

c 

1070 
C 

C 

600 
C 

C 

C 

C 

C 

C 


C 

800 


1000 

9000 


17H  TERMINAL  TIME  = .F10.3,/, 
22H  PRINTOUT  FREQUENCY  = ,15) 


IDENTIFY  VARIABLES  FOR  OUTPUT 
READ  (KIN, 1070)  LPRNTS 
READ  (KIN, 1070)  LPRNTM 
FORMAT  (6011) 


SET  THE  TIME  AND  ENTER  THE  MAIN  COMPUTATIONAL  LOOP 
TIME=T0 


START  BY  UPDATING  THE  SYSTEM 
CALL  SYSTM(TIME.X, Y,U,W,Z) 

RETURNS  X,Y  AT  T,  GIVEN  U,W,Z  OVER  (T-DEL,T) 

UPDATE  THE  MAN 

CALL  MAN ( IDUM , Y , U , Z , XMHAT , YMHAT , RES ) 

RETURNS  U OVER  (T,T+DEL)  GIVEN  Y AT  T 

CALL  INFORM  TO  DO  A PRINTOUT  IF  ONE  IS  DUE  AT  THIS  TIME 
CALL  INFORM( ISTEP , X , Y , U , XMHAT , YMHAT , RES) 

UPDATE  THE  TIME 

ISTEP=ISTEP+1 

TIME=DEL*ISTEP+TO 

IF  THE  TIME  IS  NOT  EXPIRED  DO  ANOTHER  ITERATION 
IF  (TIME+0.0001  .GE.  TEND)  GO  TO  800 
GO  TO  600 

TIME  IS  EXPIRED.  DO  OUTPUT  AND  START  AGAIN 
CALL  PRINTR(NPRINT,NPLOT,X,Z) 

GO  TO  200 

INPUT  FILE  IS  EMPTY 
WRITE  MESSAGE  AND  EXIT 
WRITE  (K0UT,9000) 

FORMAT  (IH  ,15HSTOPPING  MCARLO) 

CALL  EXIT 

END 
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MCCMP  - COMPUTATIONS  FOR  MCARLO 
INCLUDES 

1 - SUBROUTINE  SYSTM 

2 - SUBROUTINE  MAN 

3 - FUNCTION  GAUSS 


SUBROUTINE  SYSTM( T.X, Y, U. W, Z) 

INTEGRATES  SYSTEM  EQUATIONS  ONE  TIME  STEP 

HANDLES  EXTERNAL  AND  INTERNAL  SYSTEM  UPDATES 

THIS  VERSION  IS  FOR  A LINEAR  SYSTEM  USING  TRANSITION  MATRIX 


PROPAGATION 

A MORE  GENERAL  VERSION  WOULD  USE  RUNGE-KUTTA. 


1 

2 


3 


DIMENSION 

X(30).  Y(30),  U(30),  W(30),  Z(30),  ^ ^ 

DUM(26),  NST^lPdO),  ISFLAGh6),  AINT(15,15),  MDUM(16) 
XMEAN(30),  YMEAN(30),  UMEAN(30), 

XRMS(30),  YRMS(30),  URMS(30) 


INTEGER 
1 SCODE(IO) 


1 

1 

5 

8 

9 

A 

B 

C 

D 


COMMON 

/INOU/ 

/TIMES/ 

/SYSX/ 

/SYSAD/ 

/SYSY/ 

/SYSW/ 

/SYSZ/ 

/SYSINC/ 


KIN,  KOUT,  KPTR,  KPUNCH, 
KDIoK 

TIME,  DEL,  TO,  TEND,  NPRNT 

NX,  NU,  BS(15d).  AS(15,15) 
BD  15,4),  AD(l5,l5) 

NY,  DSM5,4),  CS(15,15) 

NW,  W0(4),  ES(15,4) 

NZ,  FS(15,4) 

XINC(15) 


DATA 

1 NSCODE,  LTIME,  LEND 

1 /10,  4HTIME,  AhENDS/, 

2 SCODE 

2 /1HA,  1HB,  1HC,  1HD,  1HE, 

2 1HF,  2HW0,  4HXINC,  3HINT,  5HPRINT/ 

IF  (T  .GT.  TO+0.0001)  GO  TO  100 

INITIALIZATION 

SPECIFY  THE  SYSTEM  INTERNAL  BREAKS 
READ  (KIN, 1050)  (NSTEP(I) , 1=1, NSCODE) 

1050  FORMAT  (16I5) 

WRITE  (KOUT, 1060)  ^ 

1060  FORMAT  (24H  SYSTEM  INTERNAL  BREAKS  ,17H  INDEX  CODE  NDT) 

DO  65  1=1, NSCODE 
IF  (NSTEP(I)  .LE.  0)  GO  TO  65 
WRITE  (KOUT, 1070)  I,  SCODE(I) , NSTEP(I) 

1070  FORMAT  (28X,I2,3X,A5,1X,I5) 

65  CONTINUE 

C WRITE  THE  NEXT  BATCH  OF  INPUT  CARDS  UNTIL  AN  'END'  CARD 

REWIND  KDISK 

72  READ  (KIN, 1075)  MDUM 

1075  FORMAT  (16A5) 

WRITE  (KDISK, 1075)  MDUM 
II=MDUM(1) 

IF  (lI.EQ.LEND)  GO  TO  80 
GO  TO  72 
80  TNEXT=T0 

REWIND  KDISK 
KIN1=KIN 
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C 


85 


90 

C 

100 

105 

C 

110 

120 

1130 

135 

140 

C 

150 

160 

C 

1 165 
C 

170 


1 175 

1 


C 

1 

2 

3 


IXYZ=35617 

INITIALIZE  SOME  MORE  QUANTITIES 

NW=0 

NZ=0 

DO  85  1=1,15 
XINC(I)=0.0 
DO  85  J=1,4 
DS(I,J)=0.0 
CONTINUE 

RNPTS=DEL/(TEND-TO) 

DO  90  1=1,30 

XMEAN(I)=0.0 

YMEAN(I)=0.0 

UMEAN(I)=0.0 

XRMS(I)=0.0 

YRMS(I)=0.0 

URMS(I)=0.0 

CONTINUE 

TAKE  CARE  OF  INTERNAL  SYSTEM  BREAKS 

DO  105  I=1,NSC0DE 

ISFLAG(I)=0 

IF  (NSTEP(I)  ,EQ.  0)  GO  TO  105 
ITME  = IFIX(  (T-TO+O.OOOD/DEL) 

IF  (MOD(ITME,NSTEP(I) ) .EQ.  0)  CALL  SYSNEW(I) 

CONTINUE 

TAKE  CARE  OF  EXTERNAL  SYSTEM  BREAKS 
IF  (T+0.0001  .LT.  TNEXT)  GO  TO  500 
READ  (KDISK,1130)  IDEN,  NQQ,  BRKT 
FORMAT  (A5,2X,I3,E10.0) 

IF  (IDEN. NE, LEND)  GO  TO  140 
TNEXT=1 .OE+05 
GO  TO  110 

IF  (IDEN  .NE.  LTIME)  GO  TO  150 

TNEXT=BRKT 

GO  TO  1 10 

SEARCH  THROUGH  THE  UPDATE  CODES,  SCODE(KEY) 

DO  160  KEY=1 ,NSCODE 

IF  (IDEN.EQ.SCODE(KEY))  GO  TO  170 

CONTINUE 

CODE  WAS  ILLEGAL 

WRITE  (K0UT,1165)  IDEN 

FORMAT  (23H  ILLEGAL  INPUT  CODE  OF  ,A5) 

CALL  EXIT 

DO  THE  SPECIFIED- SYSTEM  EXTERNAL  UPDATE 
ISFLAG(KEY)=1 
10=2 

KIN=KDISK 

WRITE  (K0UT.1175)  TIME,  IDEN,  NQQ 

FORMAT  (/,28H  SYSTEM  EXTERNAL  BREAK  AT  T=  ,F8 . 3 , 4X, 4HC0DE ,2X, A5 , 
4X,7HINDEX=  ,13)  . » 

GO  TO  (1,2,3,4,5,6,7,8,9,10),  KEY 

SYSTEM  DYNAMICS  - A,  B,  C,  D,  E 
NX=NQQ 

CALL  MATIO(AS,NX,NX,IO) 

GO  TO  120 
NU=NQQ 

CALL  MATIO(BS,NX,NU,IO) 

GO  TO  120 
NY=NQQ 

CALL  MATIO(CS,NY,NX,IO) 
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GO  TO  120 

4 NY=NQQ 

CALL  MATI0(DS,NY,NU,I0) 

GO  TO  120 

5 NW=NQQ 

IF  (NW.GT.O)  CALL  MATI0( ES , NX , NW, 10) 
GO  TO  120 

C DETERMINISTIC  INPUT  (F  MATRIX)  - F 

6 NZ=NQQ 

IF  (NZ.GT.O)  CALL  MATI0(FS , NX, NZ , 10) 
GO  TO  120 

C DRIVING  NOISE  - WO 

7 NW=NQQ 

IF  (NW.GT.O)  CALL  VECTI0(W0 , NW, 10) 

GO  TO  120 

C INCREMENT  TO  STATES  - XING 

8 CALL  VECTIO(XINC,NX,IO) 

GO  TO  120 

C CALL  AN  INTERNAL  UPDATE  - INT 

9 CALL  SYSNEW(NQQ) 

GO  TO  120 

C SET  PRINT  INTERVAL  - PRINT 

10  NPRNTrNQQ 
GO  TO  120 


C NO  MORE  SYSTEM  EXTERNAL  UPDATES  AT  THIS  TIME 

500  KIN=KIN1 

C COMPUTE  DISCRETE  SYSTEM  MATRICES 

IF  (ISFLAG(I)  .EQ.  0)  GO  TO  510 
CALL  DSCRT(NX,AS,DEL,AD,AINT,5) 

510  IF  (ISFLAG(1)+ISFLAG(2)  .EQ.  O)  GO  TO  520 
CALL  MMUL(AINT,BS,NX,NX,NU,BD) 

520  IF  (ISFLAG(8)  .EQ.  0)  GO  TO  540 

DO  530  1=1, NX 
X(I)=X(I)+XINC(I) 

530  CONTINUE 

C COMPUTATION  OF  TIME  AVERAGES 

540  IF  (T  .LT.  TO+0.0001)  GO  TO  630 

DO  550  1=1 , NX 

XMEAN ( I ) =XMEAN ( I ) +X ( I ) *RNPTS 
XRMS(I)=XRMS(I)+(X(I)**2)*RNPTS 
550  CONTINUE 

DO  560  1=1, NY 

YMEAN ( I ) = YMEAN ( I ) +Y ( I ) *RNPTS 
YRMS(I)=YRMS(I)+(Y(I)**2)*RNPTS 
560  CONTINUE 

DO  570  1=1, NU 

UMEAN(I)=UMEAN(I)+U(I)*RNPTS 
URMS(I)=URMS(I)+(U(I)**2)*RNPTS 
570  CONTINUE 


C 


575 

580 


COMPUTE  XRMS,  YRMS  URMS  AT  THE  END  OF  THE  RUN 
IF  (TIME+DEL+0.0001  .LT.  TEND)  GO  TO  590 
DO  575  1=1. NX 

XRMS(I)=SQRT(XRMS(I)-XMEAN(I)**2) 

CONTINUE 
DO  580  1=1, NY 

YRMS ( I ) = SORT  CYRMS ( I ) -YMEAN ( I ) **2 ) 

CONTINUE 
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DO  585  1=1, NU 

URMS(I)=SQRT(URMS(I)-UMEAN(I)»»2) 
585  CONTINUE 


C OUTPUT  THE  MEAN  AND  RMS  VALUES  OF  THE  X,  Y AND  U VECTORS 

10=3 

WRITE  (K0UT,2000) 

2000  FORMAT  ( 1?H  MEAN  OF  X VECTOR) 

CALL  VECTIO(XMEAN,NX,IO) 

WRITE  (K0UT,2010) 

2010  FORMAT  (23H  RMS  VALUES  OF  X VECTOR) 

CALL  VECTIO(XRMS,NX,IO) 

WRITE  (K0UT,2020) 

2020  FORMAT  (17H  MEAN  OF  Y VECTOR) 

CALL  VECTIO(YMEAN,NY,IO) 

WRITE  (K0UT,2030) 

2030  FORMAT  (23H  RMS  VALUES  OF  Y VECTOR) 

CALL  VECTIO(YRMS,NY,IO) 

WRITE  (K0UT,2040) 

20110  FORMAT  (17H  MEAN  OF  U VECTOR) 

CALL  VECTIO(UMEAN,NU,IO) 

WRITE  (K0UT,2050) 

2050  FORMAT  (23H  RMS  VALUES  OF  U VECTOR) 

CALL  VECTIO(URMS,NU,IO) 


C 

590 

600 


605 


INTEGRATE  SYSTEM  EQUATIONS 
DO  600  1=1, NX 
DUM(I)=0.0 
CONTINUE 

IF  (NW  .GT.  0)  CALL  VMAT2(DUM,ES,W,NX,NW,DUM) 

IF  (NZ  .GT.  0)  CALL  VMAT2( DUM, FS , Z , NX , NZ , DUM) 

DO  605  1=1, NX 

DUM(I)=DUM(I)»DEL 

CONTINUE 

CALL  VMAT2(DUM,AD,X,NX,NX,DUM) 

CALL  VMAT2(DUM,BD,U,NX,NU,X) 

CALL  VMAT1  CS,X,NY,M,DUM) 

CALL  VMAT2(DUM,DS,U,NY,NU,Y) 


C X(T+DEL)  AND  Y(T+DEL)  HAVE  JUST  BEEN  COMPUTED 

C NOW  OBTAIN  W AND  Z OVER  (T,T+DEL)  FOR  THE  NEXT  STEP 

630  IF  (T+0.0001  .GE.  TEND)  GO  TO  650 

IF  (NW  .EQ.  0)  GO  TO  640 
DO  635  1=1, NW 
C1=SQRT(W0(I)/DEL) 

W(I)  = GAUSSaXYZ)*C1 
635  CONTINUE 

640  IF  (NZ  .EQ.  0)  GO  TO  650 
DO  645  1=1, NZ 
Z(I)=FDET(I,T) 

645  CONTINUE 

650  RETURN 

END 
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SUBROUTINE  MAN ( IDUM , Y , U , Z , XH , YHAT , RES) 

INTEGRATES  MAN  EQUATIONS  ONE  TIME  STEP 

CALLS  SUBROUTINES  FOR  EXTERNAL  AND  INTERNAL  MAN  UPDATES 


1 

2 

3 


DIMENSION 


Y(30),  U(30),  Z(30),  PASTY(15,11),  P(30),  XH(30) 
VU0(4),  RES(30),  IMFLAG(32),  AVGY2(30),  AVGU(4), 
VY0(30),  FGAIN(15,15),  PASTUC(4,11) 


, YHAT(' 
AVGU2(1 


;o) 

), 


1 

1 

3 

5 

E 

F 

G 

H 

I 

J 

K 

K 

L 


COMMON 

/INOU/  KIN,  KOUT,  KPTR,  KPUNCH, 

KDISK,  IPOS,  IGOS 
/MAIN1/  NDIM,  NDIM1 , COM1(1) 

/MAIN2/  COM2(1) 

/TIMES/  TIME,  DEL  ,TO,  TEND,  NPRNT 
/MANX/  NXM,  NUM,  BM(15,4),  AM(225) 

/MANAD/  BD(15.4),  AD(15,15) 

/MANY/  NYM,  DM(15,4),  CM(225) 

/MANW/  NWM,  W0M(8),  EM(60) 

/MANZ/  NZM,  FM(60) 

/MANINC/  TD,  NPRED,  XHINC(30) 

/RATIOS/  PU(30),  VU(30),  PY(30),  VY(30),  TH(30),  ATTN(30), 
SIGMA(15, 15) 

/GAINBK/  CGN(225) 


C INITIALIZE  VECTORS  AND  MATRICES  IF  TIME  IS  LESS  THAN  TO 

IF  (TIME  .GT.  TO+0.0001)  GO  TO  100 
DO  10  1=1,30 
P(I)=0.0 
XH(I)=0.0 
YHAT(I)=0.0 
RES(I)=Y(I) 

AVGY2(I)=Y(I)*Y(I) 

10  CONTINUE 

DO  12  1=1,4 
DO  11  J=1 , 1 1 
PASTUCd,  J)=0.0 

1 1 CONTINUE 
VU0(I)=0.0 
AVGU(I)=0.0 
AVGU2(I)=0.0 

12  CONTINUE 
K0UNT=0 

CALL  IDENT(NDIM, SIGMA, 1 .OE-05) 

TC0R=1 .0 
TMEM=0.5 

ALPHA=EXP ( -DEL/TMEM) 

TPR=TMEM 
DO  13  1=1,15 
DO  13  J=1,11 
PASTYd,  J)=Y(I) 

13  CONTINUE 

C DO  EXTERNAL  AND  INTERNAL  MAN  UPDATES 

100  CALL  UPDATE(IMFLAG) 

NTOT=NXM+NUM 

IF  (IMFLAG(8)  .EQ.  1)  CALL  VADD(NXM, 1 .0,P,XHINC) 

TPR=TPR*ALPHA+DEL 

L0C=NPRED+1 

DO  110  1=1, NYM 

PASTYd, L0C)=Y(I) 

C2=PASTY( I 1 ) 

AVGY2 ( I )= ALPHA* AVGY2 (I ) +02*  02*  DEL 
C1=ABS(C2) 

C1=XGAIN(TH(I) ,0.0,01) 

C3=C1*C1*ATTN(1) 
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C1=PY(I) 

IF  (KOUNT  .LT.  NPRED)  C1=1.0E+10 
VY0(I)=(C2*C2*C1+VY(I))/C3 
RES(  I) =SQRT( VY0( I ) /DEL) *GAUSS( IDUM)+C2 
RES(I)=RES(I)-D0T3(NT0T,CM(I) ,P) 
VY0(I)=(AVGY2(I)*C1/TPR+VY(I))/C3 
110  CONTINUE 

C UPDATE  THE  KALMAN  FILTER  ESTIMATES 

CALL  MMUL ( CM , SIGMA , NYM , NTOT , NTOT , C0M1 ) 

CALL  MAT2(NYM,NT0T,C0M1 ,CM,FGAIN) 

DO  120  1=1, NYM 

FGAIN ( I . I ) =FGAIN ( 1,1) +VY0  ( I ) /DEL 
120  CONTINUE 

CALL  GMINV ( NYM , NYM , FGAIN , COM2 , MR ANK , 0 ) 

CALL  MAT5A(C0M1 , COM2, NTOT, NYM, NYM, FGAIN) 

CALL  VMAT2( P , FGAIN , RES .NTOT , NYM , P) 

CALL  MMUL ( FGAIN , CM , NTOT , NYM , NTOT , COM2 ) 

CALL  DIAG2( NTOT, COM2, COM2, -1 .0,1.0) 

CALL  MMUL (COM2, SIGMA, NTOT, NTOT, NTOT, COM 1 ) 

DO  130  1=1, NYM 
C1=VY0(I)/DEL 
DO  130  J=1,NTOT 
SIGMA(J,I)=C1*FGAIN(J,I) 

130  CONTINUE 

CALL  MAT2 ( NTOT .NYM , FGAIN , SIGMA , SIGMA ) 

CALL  MAT6S(NT0T,NT0T,C0M1 , COM2, SIGMA) 

C OBTAIN  PREDICTION  OF  CURRENT  STATE 

DO  140  1=1, NTOT 
YHAT(I)=0.0 
XH(I)=P(I) 

140  CONTINUE 

IF  (NPRED. EQ.O)  GO  TO  170 

LOC1=NPRED+1 

DO  150  L=1, NPRED 

CALL  VMAT 1( AD, XH. NTOT. NTOT, YHAT) 

CALL  VMAT2(YHAt,BD,PASTUC(1 ,L) ,NTOT,NUM,XH) 
150  CONTINUE 

170  CALL  VMAT1(CGN,XH,NUM,NXM,PASTUC(1 ,LOC)) 

DEL2=-0.5*DEL 
DO  172  1=1, NUM 
C1=ABS(VU0(I)) 

PASTUC( I , LOC) =DEL2*PASTUC(  I , LOC) 

YHAT ( I ) =DEL2* ( U ( I)+SQRT( C 1 /DEL) *GAUSS( IDUM) ) 
U(I)=U(I)+PASTUC(I,LOC) 

172  W0M(I+NWM)=C1 

II=NDIM*NXM+1 

CALL  VMAT2 ( U . CGN ( II ) , YHAT , NUM , NUM , U ) 

DO  175  1=1, NUM 

AVGU(I)=ALPHA*AVGU(I)+U(I)*DEL 

AVGU2(I)=ALPHA*AVGU2(I)+U(I)*U(I)*DEL 

C1=AVGU(I)/TPR 

VU0(I)=(AVGU2(I)/TPR-C1*C1)*PU(I)+VU(I) 

175  CONTINUE 

C PROPAGATE  SIGMA,  P 

CALL  VMAT 1( AD, P, NTOT, NTOT, YHAT) 

CALL  VMAT2( YHAt , BD , pASTUC , NTOT , NUM, P) 

CALL  VMAT 1 ( CM, XH, NYM, NTOT, YHAT) 

CALL  MMUL( AD , SIGMA , NTOT , NTOT , NTOT , COM1 ) 
NWU=NWM+NUM 
11=1 

DO  200  1=1, NWU 
C1=W0M(I)»DEL 

CALL  VSCALE(C0M2(II) ,EM(II) , NTOT, Cl) 
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II=II+NDIM 
200  CONTINUE 

CALL  MAT2 ( NTOT , NWU , COM2 , EM , SIGMA ) 

IF  (NZM.EQ.O)  GO  TO  220 
11=1 

DO  210  1=1, NZM 
C1=Z(I)»Z(I)»TC0R*DEL 
CALL  VSCALE(COM2(II) ,FM( II) , NXM, Cl ) 
II=II+NDIM 
210  CONTINUE 

CALL  MAT6S ( NXM , NZM , COM2 , FM , SIGMA ) 

220  CALL  MAT6S(NTOT, NTOT, COM1 , AD, SIGMA) 

IF  (NPRED.EQ.O)  GO  TO  235 

CALL  EQUATE(PASTY,PASTY(1 ,2) ,NYM,NPRED) 

DO  230  I=1,NPRED 

DO  230  J=1,NUM 

PASTUC( J,I)=PASTUC( J,I+1) 

230  CONTINUE 

235  KOUNT=KOUNT+1 

IF  (KOUNT.GT.NPRED)  KOUNT=NPRED 
RETURN 

C DUMMY  CALL  TO  MAT6  TO  FORCE  LOADING 

CALL  MAT6 (NTOT, NTOT, COM1 ,COM1 ,COM1) 

END 
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FUNCTION  GAUSS(DUM) 

RETURNS  A GAUSSIAN  RANDOM  VARIABLE 

WITH  ZERO  MEAN  AND  UNIT  STD  DEVIATION 

BY  SUMMING  12  UNIFORMLY  DISTRIBUTED  VARIABLES 


A=0.0 

DO  50  1=1,12 
A=A+RANF(DUM) 
50  CONTINUE 

GAUSS=A-6.0 

RETURN 

END 


60 
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MCIO  - I/O  FOR  MCARLO 
INCLUDES 

1 - SUBROUTINE  UPDATE 

2 - SUBROUTINE  INFORM 

3 - SUBROUTINE  PUTOUT 

4 - SUBROUTINE  PRINTR 


SUBROUTINE  UPDATE(IMFLAG) 

C PERFORM  EXTERNAL  UPDATES  TO  THE  MAN 


DIMENSION 

1 IMFLAG(20),  MC0DE(20),  NSTEP(32) 


1 

1 

3 

5 

8 

A 

B 

C 

E 

F 

G 

H 

I 

J 

K 

K 

L 


COMMON 

/INOU/  KIN,  KOUT,  KPTR,  KPUNCH, 

KDISK 

/MAIN1/  NDIM,  NDIM1,  COM1(1) 

/TIMES/  TIME,  DEL,  TO,  TEND,  NPRNT 
/SYSX/  NX,  NU,  B(15,4),  A(1) 

/SYSY/  NY,  D(15,4),  C(1) 

/SYSW/  NW,  W0(4),  E(1) 

/SYSZ/  NZ,  F(1) 

/MANX/  NXM,  NUM,  BM(15,4),  AM(15,15) 

/MANAD/  BD(15,4),  AD(15,15) 

/MANY/  NYM,  DM(15,4),  CM(15,15) 

/MANW/  NWM,  W0M(8),  EM(15,4) 

/MANZ/  NZM,  FM(15,4)  , ^ 

/MANINC/  TD,  NPRED,  XHINC(30)  , ^ ^ ^ 

/RATIOS/  PU(30),  VU(30),  PY(30),  VY(30) , TH(30),  ATTN(30), 
SIGMA(15,15) 

/GAINBK/  CGN(15,15) 


DATA 

1 NMCODE,  LEND,  PI,  LTIME 

1 /20,  4HENDM,  3-14159,  4HTIME/, 

2 MCODE 

2 /2HAM,  2HBM,  2HCM,  2HDM,  2HEM,  2HFM,  3HW0M,  5HXHINC, 
2 2HTD,  3HMNA,  3HMNR,  3HSNA,  3HSNR,  2HTH,  4HATTN, 

2 5HCGAIN,  5HDGAIN,  3HINT,  5HPRINT,  5HDUMMY/ 


IF  (TIME  .GT.  TO+0.0001)  GO  TO  100 


C INITIALIZATION 

C SPECIFY  THE  MAN  INTERNAL  BREAKS 

READ  (KIN, 1050)  (NSTEP(I),  1=1, NMCODE) 

1050  FORMAT  (1615) 

WRITE  (KOUT, 1060) 

1060  FORMAT  (23H  HUMAN  INTERNAL  BREAKS  ,17H  INDEX  CODE  NDT) 
DO  65  1=1, NMCODE 
IF  (NSTEP(I)  .LE.  0)  GO  TO  65 
WRITE  (KOUT, 1070)  I,  MCODE(I),  NSTEP(I) 

1070  FORMAT  (27X , 12 , 3X, A5 , IX , 15) 

65  CONTINUE 


C PARAMETER  INITIALIZATION 

TNEXT=TO 
DO  80  1=1,30 
PU(I)=0.0 
XHINC(I)=0.0 
VU(I)=0.0 
PY(I)=0.0 
VY(I)=0.0 
TH(I)=0.0 
ATTN(I)=1 .0 
80  CONTINUE 
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.85 

C 

100 

105 

C 

110 

120 

1130 

135 

140 

C 

150 

160 

C 

1165 

C 

170 
1 175 

1 

1180 

C 

1 

201 

2 

202 

3 


NPRED=0 

DO  85  I=1,NDIM 

DO  85  J=1,4 

DM(I,J)=0.0 

CONTINUE 

TD=0.0 

NWM=0 

NZMzO 


TAKE  CARE  OF  INTERNAL  MAN  BREAKS 

DO  105  I=1,NMCODE 

IMFLAG(I)=0 

IF  (NSTEP(I)  .EQ.  0)  GO  TO  105 
ITME=IFIX ( ( TIME-TO+0 . 000 1 ) /DEL) 

IF  (MOD(ITME,NSTEP(D)  .EQ.  0)  CALL  MANNEW(I) 
CONTINUE 


TAKE  CARE  OF  EXTERNAL  MAN  BREAKS 
IF  (TIME+0.0001  .LT.  TNEXT)  GO  TO  500 
READ  (KIN, 1130)  IDEN,  NQQ,  BRKT 
FORMAT  (A5,2X,I3,E10.0) 

IF  (IDEN. NE. LEND)  GO  TO  140 
TNEXT=1 .OE+05 
GO  TO  1 10 


IF  (IDEN.NE.LTIME)  GO  TO  150 

TNEXT=BRKT 

GO  TO  110 


SEARCH  THROUGH  THE  UPDATE  CODES 
DO  160  KEY=1,NMCODE 
IF  ( IDEN. EQ.MC0DE( KEY))  GO  TO  170 
CONTINUE 


CODE  WAS  ILLEGAL 

WRITE  (K0UT,1165)  IDEN 

FORMAT  (23H  ILLEGAL  INPUT  CODE  OF  ,A5) 

CALL  EXIT 


DO  THE  SPECIFIED  MAN  EXTERNAL  UPDATE 
IMFLAG(KEY)=1 

WRITE  (K0UT,1175)  TIME,  IDEN,  NQQ 

FORMAT  (/,27H  HUMAN  EXTERNAL  BREAK  AT  T=  ,F8 .3 , 4X , 4HCODE 
4X,7HINDEX=  ,13) 

IF  (NQQ.LE.O  .AND.  KEY.LE.7)  WRITE  (KOUT,1180) 

FORMAT  (21H  HUMAN  MODEL  = SYSTEM) 

10=2 

GO  TO  (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19), 

SYSTEM  DYNAMICS  --AM,  BM,  CM,  DM,  EM 
NXM=NQQ 

IF  (NXM.LE.O)  GO  TO  201 
CALL  MATI0(AM,NXM,NXM,I0) 

GO  TO  120 
NXMzNX 

CALL  EQUATE(AM,A,NXM,NXM) 

GO  TO  120 
NUM=NQQ 

IF  (NUM.LE.O)  GO  TO  202 
CALL  MATIO(BM,NXM,NUM,IO) 

GO  TO  120 
NUMrNU 

CALL  EQUATE(BM,B,NXM,NUM) 

GO  TO  120 
NYMrNQQ 

IF  (NYM.LE.O)  GO  TO  203 
CALL  MATIO(CM,NYM,NXM,IO) 


2X,A5, 

KEY 
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GO  TO  120 

203  NYM=NY 

CALL  EQUATE(CM,C,NYM,NXM) 

GO  TO  120 

4 NYMrNQQ 

IF  (NYM.LE.O)  GO  TO  204 

CALL  MATIO(DM,NYM,NUM,IO) 

GO  TO  120 

204  NYM=NY 

CALL  EQUATE(DM,D,NYM,NXM) 

GO  TO  120 

5 NWM=NQQ 

IF  (NWM.LE.O)  GO  TO  205 

CALL  MATIO(EM,NXM,NWM,IO) 

GO  TO  120 

205  NWM=NW 

IF  (NW.GT.O)  CALL  EQUATE ( EM, E,NXM,NWM) 
GO  TO  120 

C DETERMINISTIC  INPUT  (FM  MATRIX)  - FM 

6 NZMzNQQ 

IF  (NZM.LE.O)  GO  TO  206 

CALL  MATIO(FM,NXM,NZM,IO) 

GO  TO  120 

206  NZMzNZ 

IF  (NZ.GT.O)  CALL  EQUATE(FM,F,NXM,NZM) 
GO  TO  120 

C DRIVING  NOISE  - WOM 

7 NWMzNQQ 

IF  (NWM.LE.O)  GO  TO  207 
CALL  VECTIO(WOM,NWM,IO) 

GO  TO  120 

207  NWM=NW 

IF  (NW.GT.O)  CALL  EQUATE(WOM,WO,NWM, 1 ) 
GO  TO  120 

C INCREMENT  TO  STATE  ESTIMATES  - XHINC 

8 CALL  VECTIO(XHINC,NXM,IO) 

GO  TO  120 

C TIME  DELAY  - TD 

9 CALL  VECTIO(TD,1,IO) 

NPRED=IFIX( TD/DEL+0 .5001) 

TD=DEL*NPRED 

GO  TO  120 

C NOISES  - MNA,  MNR,  SNA,  SNR 

10  CALL  VECTIO(VU,NUM,IO) 

GO  TO  120 

11  CALL  VECTIO(PU,NUM,IO) 

DO  211  1=1, NUM 
PU(I)=PI*10.0**(PU(I)/10.0) 

211  CONTINUE 

GO  TO  120 

12  CALL  VECTIO(VY,NYM,IO) 

GO  TO  120 

13  CALL  VECTIO(PY,NYM,IO) 

DO  213  1=1, NYM 
PY(I)=PI*10.0**(PY(I)/10.0) 

213  CONTINUE 

GO  TO  120 

C THRESHOLDS  AND  ATTENTION  - TH,  ATTN 

14  CALL  VECTIO(TH,NYM;IO) 

GO  TO  120 
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15  CALL  VECTIO(ATTN,NYM,IO) 

GO  TO  120 

C CONTINUOUS  AND  DISCRETE  CONTROL  GAINS  - CGAIN,  DGAIN 

16  CALL  MATIO(CGN,NUM,NUM+NXM,IO) 

GO  TO  120 

17  CALL  MATIO(CGN,NUM,NUM+NXM,IO) 

GO  TO  120 

C CALL  AN  INTERNAL  MAN  UPDATE  - INT 

18  CALL  MANNEW(NQQ) 

GO  TO  120 

C PRINT  INTERVAL  - PRINT 

19  NPRNT=NQQ 
GO  TO  120 

C NO  MORE  MAN  UPDATES  AT  THIS  TIME 

500  CONTINUE 

C UPDATE  VARIOUS  MAN  PARAMETERS 

NTOT=NXM+NUM 

C COMPUTE  DISCRETE  MAN  MATRICES 

IF  (IMFLAG(1)+IMFLAG(2)  .EQ.  0)  GO  TO  520 

DO  510  1=1 ,NUM 

II=NXM+I 

DO  505  J=1,NXM 

AD(II, J)p0.0 

505  CONTINUE 

DO  506  K=1 ,NUM 
BD(II,K)=0.0 

506  CONTINUE 
BD(II,I)=1 .0 

510  CONTINUE 

CALL  DSCRT(NXM,AM,DEL,AD,COM1 ,5) 

CALL  MMUL(COM1 ,BM,NXM,NXM,NUM,BD) 

CALL  EQUATE ( AM ( 1 , NXM+1 ) ,BM,NXM,NUM) 

C INCORPORATE  NEW  CGAINS 

520  IF  (IMFLAG(16) .EQ.O)  GO  TO  5^0 

CALL  MSCALE( AM (NXM+1 , 1 ) . CGN . NUM, NTOT , -1 . 0) 

CALL  DSCRT(NTOT,AM,DEL.CGN,COM1 ,5) 

CALL  MMUL(AM(NXM+1 ,1) , COM 1 , NUM, NTOT, NTOT, CGN) 

DO  525  1=1, NUM 
II=I+NXM 
DO  525  J=1iNTOT 
AM(II,J)=0.0 
CGN(I, J)=-CGN(I, J)/DEL 
525  CONTINUE 

WRITE  (K0UT,5250) 

5250  FORMAT  (37H  EQUIVALENT  DISCRETE  GAINS  GENERATED  ) 

10=3 

CALL  MATIO(CGN,NUM,NTOT, 10) 

C UPDATE  EM  AND  CM 

5^10  CALL  MMUL(BD,CGN(1, NXM+1), NTOT, NUM, NUM, EM(1,NWM+D) 

CALL  EQUATE(CM( 1 , NXM+1 ) , DM, NYM, NUM) 

DO  550  J=1 ,NUM 

JJ=J+NXM 

JQ=J+NWM 

DO  550  1=1, NTOT 

EM(I, JQ)=0.5*EM(I,JQ) 

ADa,JJ)=BD(I,J)-DEL*EM(I,JQ) 

550  CONTINUE 

IF  (IMFLAG(5)+IMFLAG(7) .EQ.O)  GO  TO  570 
IQ=NXM+1 

IF  (NWM.LE.O)  GO  TO  570 
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DO  560  J=1,NWM 
DO  560  I=IQ,NT0T 
560  EM(I,J)=0.0 

570  CONTINUE 

RETURN 

END 
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SUBROUTINE  INF0RM( ISTEP , X , Y , U . XH , YH , RES ) 
C DO  OUTPUT  FOR  A SINGLE  TIME  STEP 

DIMENSION 

1 X(1),  Y(1),  U(1),  XH(1),  YH(1),  RES(I), 

COMMON 

1 /INOU/  KIN,  KOUT.  KPTR.  KPUNCH, 

1 KDISK,  IPOS.  IGOS 

5 /TIMES/  TIME,  DEL,  TO,  TEND,  NPRNT 

6 /INFO/  NREC,  NPRINT,  NPLOT,  LPP(20,6), 

6 SMIN(21),  SMAXC21) 

C CHECK  IF  TIME  FOR  SOME  OUTPUT 

IF  (MOD( ISTEP, NPRNT)  .NE.  0)  RETURN 

C DO  SOME  INITIALIZATION 

IF  (TIME  .GT.  TO+0.0001)  GO  TO  100 
NPRINT=1 
NPLOT=1 
REWIND  IPOS 
REWIND  IGOS 
DO  8 1=1,21 
SMIN(I)=1 .OE+20 
SMAX(I)=-1 .OE+20 
8 CONTINUE 

DO  10  1=1,6 
IP(I)=0 
IGU)  = 0 
DO  9 J=1,20 
L1=LPP(J,I) 


IF 

(LI .EQ.1 

.OR. 

LI .EQ.2) 

IP(] 

[)=ip(: 

[)  + 1 

IF 

(LI .EQ.2 

.OR. 

L1.EQ.3) 

IG(] 

[)  = IG(j 

[)  + 1 

9 CONTINUE 
NPRINT=NPRINT+IP(I) 
NPLOT=NPLOT+IG(I) 

10  CONTINUE 
NREC=0 

IF  ( ISTEP. EQ.O)  RETURN 

C DO  OUTPUT  FOR  THE  CURRENT  TIME 

100  DUM1(1)=TIME 

DUM2(1)=TIME 

IF  (SMIN(1) .GE.TIME)  SMIN(1)=TIME 

SMAX(1)=TIME 

LOC1=1 

L0C2=1 

DO  160  J=1,6 
DO  160  1=1,20 
L1=LPP(I, J) 

IF  (LI .GE.3)  GO  TO  1 40 

IF  (L1.LE.0)  GO  TO  160 

LOC1=LOC1+1 

IF  (J.EQ.1)  DUM1(LOC1)=X(I) 

IF  (J.EQ.2)  DUM1(LOC1)=Y(I) 

IF  (J.EQ.3)  DUM1(L0C1)=U(I) 

IF  (J.EQ.4)  DUM1(L0C1)=XH(I) 

IF  (J.EQ.5)  DUM1(L0C1)=YH(I) 

IF  (J.EQ.6)  DUM1(LOC1)=RES(I) 

140  IF  (LI .EQ.1)  GO  TO  160 

L0C2=L0C2+1 

IF  (J.EQ.1)  DUM2(L0C2)=X(I) 

IF  (J.EQ.2)  DUM2(L0C2)=Y(I) 

IF  (J.EQ.3)  DUM2(L0C2)=U(I) 

IF  (J.EQ.4)  DUM2(L0C2)=XH(I) 

IF  (J.EQ.5)  DUM2(L0C2)=YH(I) 


DUMK21),  DUM2(21) 


IP(6),  IG(6), 
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IF  (J.EQ.6)  DUM2(L0C2)=RES(I) 

150  C1=DUM2(L0C2) 

IF  (SMIN(L0C2) .GE.C1)  SMIN( L0C2) =C1 
IF  (SMAX(L0C2) .LE.C1)  SMAX( L0C2 ) =C1 
160  CONTINUE 

IF  (L0C1.GT.1)  CALL  PUT0UT(DUM1 ,NPRINT, IPOS) 
IF  (L0C2.GT.1)  CALL  PUT0UT(DUM2 , NPLOT, IGOS) 
NREC=NREC+1 
RETURN 

END 
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SUBROUTINE  PUTOUT( DUM , NVAR , IDISK) 

C SUBROUTINE  TO  SAVE  OUTPUT  ON  A FILE 

DIMENSION 
1 DUM(1) 

WRITE  (IDISK)  (DUM(I),  1=1, NVAR) 
RETURN 

END 
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SUBROUTINE  PRINTR(NPRINT,NPL0T,DUM1 ,DUM2) 
C PRINT  THE  OUTPUT  AT  THE  END  OF  A RUN 


1 


DIMENSION  . ^ 

DUMId),  DUM2(1),  GRAPH(1350),  TITLE(6),  LET(11) 


1 

1 

2 

2 

3 

6 

6 


COMMON 

/INOU/ 

/PL0T1/ 

/MAIN1/ 

/MAIN2/ 

/INFO/ 


KIN,  KOUT,  KPTR.  KPUNCH, 

KDISK,  IPOS,  IGOS 

NV,  nA,  NCP1&,  LW,  XL,  XH,  YL , YH,  NXES,  NDIR, 
NGLV,  NGLH,  BSYM,  GSYM,  PSYM,  ND1 , ND2,  NOUT 
NDIM^^NDIMI,  STORE(I) 

NREC,  II,  12,  LPP(20,6),  IP(6),  IG(6), 
SMIN(21),  SMAX(21) 


1ST, 


DATA 
1 TITLE 

1 /8H  STATE  , 8H  OUTPUT  , 8HCONTROL  , 8H  XMHAT  , 
1 8H  YMHAT  , 8H  INOVAT  / 

IF  (NPRINT.EQ.1)  GO  TO  51 
REWIND  IPOS 
DO  10  1=1, NREC 

READ  (IPOS)  (DUMI(KK),  KK=1,NPRINT) 

II=I 

DO  9 L=1,NPRINT 
STORE(II)=DUM1(L) 

9 II=II+NREC 

10  CONTINUE 


30 


1035 


1045 


40 

50 


IBEG=NREC+1 
DO  50  1=1,6 
M=0 


DO  30  L=1,20 
IQ=LPP(L,I) 

IF  (IQ.EQ.O  .OR.  IQ.EQ.3)  GO  TO  30 

M=M+1 

LET(M)=L 

CONTINUE 

IF  (M.EQ.O)  GO  TO  50 
CALL  PAGEFD(KOUT, 1) 

WRITE  (KOUT, 1035)  (TITLE(I),  LET( J) , J=1 ,M) 

FORMAT  (1H  ,3X,4HTIME,2X,10(a8,I2,2X)) 

LIM1=IBEG 

LIM2=IBEG+(M-1)*NREC 
DO  40  L=1 ,NREC 

WRITE  (KOUT, 1045)  STORE(L) , (STORE(J) , J=LIM1 , LIM2 , NREC) 
FORMAT  (1H  ,F7.2,10(2X,1PE10.3)) 

LIM1=LIM1+1 

LIM2=LIM2+1 

CONTINUE 

IBEG=IBEG+M*NREC 


51  IF  (NPLOT.EQ.1)  RETURN 

REWIND  IGOS 


DO  60  1=1. NREC 

READ  (IGOS)  (DUM2(KK) ,KK=1 ,NPLOT) 
II=I 

DO  59  L=1,NPLOT 
STORE(II)=DUM2(L) 

59  II=II+NREC 

60  CONTINUE 
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ND2=NREC 
XH=SMAX( 1 ) 

XL=SMIN(1) 

M=1 

DO  70  1=1,6 
DO  70  L=1,20 
IQ=LPP(L,I) 

IF  (IQ.LT.2)  GO  TO  70 
M=M+1 

YH=SMAX(M) 

YL=SMIN(M) 

IBEG=(M-1)»NREC+1 

CALL  KPLOT(GRAPH,STORE,STORE(IBEG) ,0,0,0, 0,0) 
WRITE  (KOUT, 1082)  TITLE(I),  L 
1082  FORMAT  (/,1H  ,A8,I3) 

70  CONTINUE 

RETURN 

END 
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SUBROUTINE  KPLOT ( W , X , Y , NTAPE , IX , lY , NVAR , Y 1 ) 


COMMON  /PL0T1/ 

1 NV , NH , NCPW , LW , VLH ( 4 ) , NXES , NDIR , 1ST , NGLV , NGLH , BSYM , GSYM , 

2 PSYM , NDIM 1 , NDIM2 , NO 

DIMENSION  W(1)  ,X(1)  ,Y(1)  .YKNVAR)  ,STORE(70)  ,Q(4)  ,IPX(4)  ,K(3) 


EQUIVALENCE 

EQUIVALENCE 


(Q(1)  ,XL1K(Q(2)  ,XH1)  ,(Q(3)  ,YL1)  ,(Q(4)  ,YH1) 
USC,K(D)  ,(JSC,K(3)) 


DATA  IPX/3, 4, 1 ,2/ 

IF  (NH.GT.121)  NH=121 

C NCPW  IS  THE  NUMBER  OF  CHARACTERS  PER  WORD 

C (60  BIT  WORD  6 BIT  DISPLAY  CODE  ON  CDC) 

NCPWrIO 
LW=NH/NCPW+1 

IF  ( (IST/10) .GT.O)  NCOUNTrO 
NCOUNT=NCOUNT+1 
IF  (NCOUNT.EQ. 10)  IST=1 
L=1 

DO  10  1=1,4 

Q(I)=-1 .0E08*(-1)**I 

K(L)=1 

IF  (VLH(L) .EQ.VLH(L+1 ) ) GO  TO  10 

K(L)=0 

Q(I)=VLH(I) 

10  IF  I.EQ.2)  L=3 

IF  ( NTAPE. EQ.O)  GO  TO  1200 


C SKIP  THIS  PART  IF  PLOTTING  FROM  CORE 

IFLAG=0 
GO  TO  40 
1600  IFLAG=1 

40  NN=0 

REWIND  NTAPE 
50  READ  (NTAPE)  Y1 

C GO  TO  2800  ON  EOF 

IF  (EOF(NTAPE))2800, 100 
100  NN=NN+1 

IF  (NN.LT.NDIM1)  GO  TO  50 
IF  (IFLAG.EQ.1)  GO  TO  1700 
IF  (ISC+ JSC. EQ.O)  GO  TO  1710 
300  NN=NN+1 

IF  (ISC. EQ.O)  GO  TO  600 
XL1=AMIN1(XL1 ,Y1(IX)) 
XH1=AMAX1(XH1,Y1(IX)) 

IF  (JSC. EQ.O)  GO  TO  200 
600  YL1=AMIN1(YL1 ,Y1(IY)) 

YH1=AMAX1(YH1 ,Y1(IY)) 

200  READ  (NTAPE)  Y1 

IF  (EOF(NTAPE) ) 1600,210 
210  IF  (NN-NDIM2)  300,300,1600 


C RESUME  HERE 

1200  IF  (ISC. EQ.O)  GO  TO  1400 

DO  1300  I=NDIM1 ,NDIM2 
XL1=AMIN1(XL1  ,X(D) 

1300  XH1=AMAX1(XH1  ,X(D) 

1400  IF  (JSC. EQ.O)  GO  TO  1700 

DO  1500  I=NDIM1 ,NDIM2 
YL1=AMIN1(YL1,Y(I)) 

1500  YH1=AMAX1(YH1,Y(D) 

1700  IF  (ISC.EQ.1)  CALL  ADJUST(XH1 , XL1 ) 

IF  (JSC.EQ.1)  CALL  ADJUST( YH1 , YL1 ) 
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1710 

1720 


1740 

1760 


1780 


1800 


2200 

C 


2400 


C 

2500 

2600 

C 

2700 


C 

2800 

2900 


IF  (NDIR/10)  1720,1740,1720 

TMPrXLI 

XL1=XH1 

XH1=TMP 

IF  (NDIR-10*(NDIR/10))  1760,1780,1760 

TMP=YL1 

YL1=YH1 

YH1=TMP 

J=7*(NC0UNT-1)+1 
IF  (J.EQ.1)  CALL  QINIT(W) 

STORE(J)=PSYM 
DO  1800  1=1,4 
IF  (NXES.EQ.O)  L=I+J 
IF  (NXES.GT.O)  L=IPX(I)+J 
ST0RE(L)=Q(I) 

ST0RE( J+5 ) = ( NH- 1 ) / (ST0RE( J+2 ) -ST0RE( J+1 ) ) 

ST0RE(J+6)  = (NV-1)/(ST0REU+4)-ST0REU+3)) 

IF  (NTAPE.EQ.O)  GO  TO  2500 

SKIP  THIS  PART  IF  PLOTTING  FROM  CORE 
DO  2400  I=NDIM1 ,NDIM2 

IF  (NXES.EQ.O)  CALL  KPLOTC( STORE( J) , W, Y1 ( IX) , Y1 ( lY) ) 

IF  (NXES.GT.O)  CALL  KPLOTC(STORE( J) ,W,Y1 (lY) ,Y1 (IX) ) 

READ  (NTAPE)  Y1 
IF  (EOF(NTAPE))  2800,2700 

SKIP  THIS  PART  IF  PLOTTING  FROM  A FILE 
DO  2600  I=NDIM1 ,NDIM2 

IF  (NXES.EQ.O)  CALL  KPLOTC( STORE( J) , W , X( I) , Y( I) ) 

IF  (NXES.GT.O)  CALL  KPLOTC(sTORE(j) ,W,y(i) ,x(l)) 

RESUME  HERE 

IF  (dST-10*(IST/10))  .GT.O)  CALL  QPRINT(W, NO, NCOUNT, STORE) 
RETURN 

ERROR  MESSAGE 
WRITE  (NO, 2900) 

FORMAT (//32H  INSUFFICIENT  DATA  ON  INPUT  FILE,/, 

1H  ,28H  PLOTTING  ROUTINE  TERMINATED) 

RETURN 

END 
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SUBROUTINE  ADJUST(XH1 ,XL1 ) 

IF  (XH1.EQ.XL1)  XL1=0.9*XL1-10.0 
A=IFIX(100.0+ALOG10(XH1-XL1))-100.0 
XH1T=XH1»10.0*»(1 .0-A) 
XL1T=XL1»10.0»»( 1 .0-A) 

IF  (XH1T  .GE.  0.0)  XH1T=IFIX(XH1T+0. 
XH1T=IFIX(XH1T) 

IF  (XL1T  .LE.  0.0)  XL1T=IFIX(XL1T-0. 
XL1T=IFIX(XL1T) 

XH1=XH1T*10.0»»(A-1 .0) 
XL1=XL1T»10.0»*(A-1 .0) 


9) 

9) 


RETURN 

END 
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SUBROUTINE  QINIT( IMAGE) 

COMMON  /PL0T1/ 

1 NV , NH , NC  PW , LW , Q ( 4 ) , NXES , NDIR , 1ST , NGLV , NGLH , BSYM , GSYM 

DIMENSION  IMAGE(I) 

DATA  IBLNK/10H 

N=LW*NV 
DO  100  1=1, N 

100  IMAGE(I)=IBLNK 
DO  101  1=1, NH 

CALL  QPL0T( IMAGE ,1,1, BSYM) 

101  CALL  QPL0T( IMAGE, I, NV, BSYM) 

DO  102  1=1, NV 

CALL  QPL0T( IMAGE, 1,1, BSYM) 

102  CALL  QPLOTC IMAGE, NH, I, BSYM) 

1800  IF  (NGLV.EQ.O)  GO  TO  2000 

NGLV1=NGLV+1 

NH1=NH-1 

DO  1900  I=NGLV1 ,NH1 ,NGLV 
DO  1900  J=1,NV 

1900  CALL  QPLOT(IMAGE,I.J,GSYM) 

2000  IF  (NGLH.EQ.O)  RETURN 

NGLH1=NGLH+1 
NV1=NV-1 

DO  2100  I=NGLH1 ,NV1 ,NGLH 
DO  2100  J=1.NH 

2100  CALL  QPLOTCIMAGE, J,I,GSYM) 

RETURN 

END 
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SUBROUTINE  KPL0TC(W, IMAGE, 

DIMENSION  W(1),  IMAGE (1) 
COMMON  /PL0T1/  NV,  NH 

J=(X-W(2) )*W(6)+1 .5 

IF  (( J.LE.O) .OR. (J.GT.NH)) 

I=NV-IFIX( (Y-W(4))*W(7)+0. 

IF  ((I.LE.O).OR.(I.GT.NV)) 

CALL  QPLOT  (IMAGE , J , I , W(1 ) 

RETURN 

END 


,Y) 


RETURN 

) 

RETURN 
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SUBROUTINE  QPLOT( IMAGE , J , I ,SYM) 

COMMON  /PL0T1/  NV,  NH,  NCPW,  LW 
DIMENSION  IMAGE(I) 

II=J/NCPW 

L=J-NCPW*II 

11=11+1 

IF  (L)  101,101,102 

101  L=NCPW 
11=11-1 

102  IW=II+(I-1)*LW 

CALL  PLACEdMAGE(IW)  ,L,SYM,  1) 

103  RETURN 
END 
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SUBROUTINE  PLACE(A,N,B,M) 

THE  MTH  CHAR  OP  B REPLACES 
THE  NTH  CHARACTER  OF  A 

CHAR  POSITIONS  ARE  1 TO  10  FROM  LEFT  TO  RIGHT 

COMMON/INOU/KIN , KOUT 
INTEGER  A,  B,  BX,  BY 
DATA  MASK/7 78/ 

CHECK  FOR  VALID  ARGUMENTS 
IF  (N.GT.10  .OR.  M.GT.10)  GO  TO  900 
IF  (N.LT.1  .OR.  M.LT.1)  GO  TO  900 

NULL  ALL  BUT  THE  MTH  CHAR  OF  B,  PUT  IT  IN  BX 
NULL  THE  NTH  CHAR  OF  A 
NSHFT=60-6*N 
MSHFT=60-6*M 

MASKBY  = SHIFT ( MASK, NSHFT) 

MASKB  = SHIFT ( MASK, MSHFT) 

MASKA  = COMPL (MASKBY) 

A = AND (A, MASKA) 

BX  = AND (B, MASKB) 

SHIFT  THE  MTH  CHAR  OF  BX  TO  THE  NTH  POSITION 
PUT  IT  IN  BY  AND  NULL  ALL  BUT  THE  NTH  CHAR 
MNSHFT=6*(M-N) 

BY  = SHIFT (BX,MNSHFT) 

BY  = AND (BY, MASKBY) 

C COMBINE  A AND  BY 

A = OR(A,BY) 

RETURN 

C N OR  M OUT  OF  BOUNDS 

900  WRITE  (KOUT, 1900) 

1900  FORMAT  (1H  .20HERROR  IN  SUBR.  PLACE,/, 

1 2i|H  N OR  M IS  OUT  OF  BOUNDS) 

CALL  EXIT 

END 
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SUBROUTINE  QPRINT(IMAGE , NO , NCOUNT .STORE) 


DIMENSION  IMAGE(I),  STORE(I) 
COMMON  /PL0T1/  NV,  NH,  NCPW,  LW 


1 10 


120 


130 

150 

102 

103 


CALL  PAGEFD(N0,1) 

DO  110  1=1, NCOUNT 
IB=7*(I-1)+1 

WRITE  (NO, 102)  STORE(IB+1 ) ,STORE(IB) ,STORE(IB) ,STORE(IB+2) 

NCANT=NV-NCOUNT 

IA=1 


DO  150  1=1 ,NV 
IB=I*LW 

IF  ( I. GT. NCOUNT)  GO  TO  120 
IBASE=(I-1)*7+1 

WRITE  (NO, 103)  STORE(IBASE) ,STORE(IBASE+4) ,(IMAGE(J) ,J=IA,IB) 
GO  TO  150 

IF  (I.GT.NCANT)  GO  TO  13O 
WRITE  (NO, 105)  (IMAGE(J) ,J=IA,IB) 

GO  TO  150 

IBASE=(I-1-NCANT)*7+1 

WRITE  (no, 103)  STORE(IBASE) ,STORE(IBASE+3) ,(IMAGE(J) ,J=IA,IB) 
IA=IA+LW 


FORMAT ( 
FORMAT ( 
FORMAT ( 
RETURN 
END 


1H  ,11X,1PE10.3,1X,A1,77X,A1 
1H  ,A1, 1PE9.2, IX, 12A10,A1) 

1H  ,11X,12A10,A1) 


1PE10. 


3) 
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C 


10 

C 


20 

30 


50 

60 

70 

80 


SUBROUTINE  DSCRT( N , A , DEL , EA , EAINT , NT) 

DIMENSION  A(1) ,EA( 1 ) ,EAINT( 1 ) ,COEF(30) 

SETS  EA=EXP(A*DEL) ,EAINT=INTEGRAL  EA  0 TO  DEL 
COMMON/MAIN1 /NDIM , NDIM1 
NN=N*NDIM 
NTM1=NT-1 
C0EF(NT)=1 . 

DO  10  I=1,NTM1 
II=NT-I 

C0EF(II)=DEL*C0EF(II+1)/FL0AT(I) 

NT  MUST  BE  AT  LEAST  3 

11=1 

DO  30  1=1, N 
DO  20  J=I,NN,NDIM 
EAINT(J)=A(J)*COEF(D 
EAINT( II ) =EAINT( II )+C0EF( 2 ) 

II=II+NDIM1 
DO  60  L=3,NT 
T1=C0EF(L) 

CALL  MMUL(A,EAINT.N,N,N,EA) 

IF(L.EQ.NT)GO  TO  70 
11=1 

DO  60  1=1 ,N 
DO  50  J=I,NN,NDIM 
EAINT(J)=EA( J) 

EAINT(II)=EAINT(II)+T1 

II=II+NDIM1 

DO  80  11=1 ,NN,NDIM1 

EA(II)=EA(II)+T1 

CONTINUE 

RETURN 

END 
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SUBROUTINE  GMINV ( NR , NC , A , U , MR , MT ) 
DIMENSION  A(1) ,U(1) ,S(30) 

COMMON/MAIN 1/  NDIM,NDIM1 

COMMON/INOU/  KIN.KOUT 

T0L=1 .E-12 

MR=NC 

NRM1=NR-1 

TOL1=1 .E-20 

JJ=1 

DO  100  J=1,NC 
FAC=DOT(NR,A( JJ) ,A(JJ)) 

JM1=J-1 
JRM=JJ+NRM1 
JCM=JJ+JM1 
DO  20  I=JJ,JCM 
20  U(I)=0. 

U(JCM)=1 .0 
IF(J.EQ.I)  GO  TO  54 
KK=  1 

DO  30  K=1.JM1 

IF(S(K) .EQ. 1 .0)  GO  TO  30 

TEMP=-DOT(NR,A(JJ) ,A(KK)) 

CALL  VADD(K,TEMP,U(JJ) ,U(KK)) 

30  KK=KK+NDIM 

DO  50  L=1,2 
KK=  1 

DO'50  K=1,JM1  ■ 

IF(S(K) .EQ.O. ) GO  TO  50 
TEMP=-DOT(NR,A(JJ) ,A(KK)) 

CALL  VADD(NR,TEMP,A(JJ) ,A(KK)) 

CALL  VADD(K,TEMP,U(JJ) ,U(KK)) 

50  KK=KK+NDIM 

TOL1=TOL*FAC 
FAC=DOT(NR,A(JJ) ,A(JJ)) 

54  IF(FAC.GT.TOLI)  GO  TO  70 

DO  55  I=JJ,JRM 

55  A(I)=0. 

S(J)=0. 

KK=  1 

DO  65  K=1,JM1 

IF(S(K) .EQ.O. ) GO  TO  65 

TEMP=-DOT(K,U(KK) ,U(JJ)) 

CALL  VADD(NR,TEMP,A(JJ) ,A(KK)) 

65  KK=KK+NDIM 

FAC=DOT(J,U(JJ) ,U(JJ) ) 

MR=MR-1 
GO  TO  75 
70  S(J)=1.0 

KK=  1 

DO'72  K=1,JM1 

IF(S(K) .EQ. 1 . ) GO  TO  72 

TEMP=-DOT(NR,A( JJ) ,A(KK) ) 

CALL  VADD(K,TEMP,U(JJ) ,U(KK)) 

72  KK=KK+NDIM 

75  FAC=1 ./SQRT(FAC) 

DO  80  I=JJ,JRM 
80  A(I)=A(I)*FAC 

DO  85  I=JJ,JCM 
85  U(I)=U(I)*FAC 

100  JJ=JJ+NDIM 

IFCMR.EQ. NR.OR.MR.EQ.no  GO  TO  120 
IF(MT.NE.O)  WRITE  (K0UT,110)  NR,NC,MR 
110  F0RMAT(I3, 1HX,I2,8H  M RANK, 12) 

120  NEND=NC*NDIM 

JJ=1 
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125 


DO  135  J=1,NC 
DO  125  1=1, NR 
II=I-J 
S(I)=0. 

DO  125  KK=JJ,NEND,NDIM 
S(I)=S(I)+A(II+KK)*U(KK) 
II  = J 

DO  130  1=1, NR 
U(II)=S(I) 

II=II+NDIM 

JJ=JJ+NDIM1 

RETURN 

END 
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SUBROUTINE  MAT2(N1 ,N2,X,Y,Z) 
C ZrXY  X,Y=N1»N2,Z=Z' 

C Z AND  Y CAN  BE  EQUIVALENT 

DIMENSION  X(1) ,Y(1) ,Z(1) 
C0MM0N/MAIN1/NDIM,NDIM1 
NN2=N2*NDIM 
11=1 

DO  10  1=1, N1 
IJ=II 

DO  5 J=I.N1 

Z(IJ)=D0T2(NN2,X(I) ,Y(J)) 

5 IJ=IJ+NDIM 

J=II 
IJ=J 

3 IJ=IJ-NDIM 

IF(IJ.LT.I)  GO  TO  10 
J=J-1 

Z(IJ)=Z(J) 

GO  TO  3 

10  II=II+NDIM1 

RETURN 
END 
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SUBROUTINE  MAT5A(X, Y,N1 ,N2,N3,Z) 

C Z=XT»Y  X=N2»N1,  Y=N2»N3 

DIMENSION  X(1) ,Y(1) ,Z(1) 

C0MM0N/MAIN1/NDIM 
N1M1=N1-1 
NN3=N3*NDIM 
DO  1 I=1,NN3,NDIM 
II=I+N1M1 
DO  1 J=I,II 
1 Z(J)=0.0 

ENTRY  MAT5AS 
NN3=N3*NDIM 
DO  10  K=1,N2 
KK=K 

DO  8 1=1, N1 
C 1 ”X( KK) 

IF(CI.NE.O.O)  CALL  VADD1 (NN3 , Cl , Z( I) , Y(K) ) 
8 KK=KK+NDIM 

10  CONTINUE 

RETURN 
END 
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SUBROUTINE  MAT6(N1 ,N2,X, Y,Z) 

C Z=X*Y  , WHERE  X=N1*N2,  Y=N1*N2,  Z=Z'=N1*N1 

DIMENSION  X(1)  , Y(1)  , Z(1) 

COMMON  /MAIN1/  NDIM,  NDIM1 

NN1=N1*NDIM 
DO  1 1=1 ,N1 
DO  1 J=I,NN1,NDIM 
Z(J)=0.0 
1 CONTINUE 

ENTRY  MAT6S 
C Z=Z+X*Y' 

NN2=N2*NDIM 

NN1=N1*NDIM 

DO  6 K=1,NN2,NDIM 

KK=K-1 

J=1 

DO  6 1=1, N1 
C1=Y(I+KK) 

IF  (C1.NE.0.0)  CALL  VADD(I,C1 ,Z(J) ,X(K) ) 
J=J+NDIM 
6 CONTINUE 

IF  (N1.EQ.1)  RETURN 

NN2=NDIM1+1 

DO  10  K=NN2,NN1 ,NDIM1 

I = K 

J=K 

8 1=1-1 

J=J-NDIM 

Z(J)=Z(I) 

IF  (J.GT.NDIM)  GO  TO  8 
10  CONTINUE 

RETURN 
END 


84 


Report  No.  3^63 


Bolt  Beranek  and  Newman  Inc. 


SUBROUTINE  MMUL(X,Y,N1,N2,N3,Z) 

DIMENSION  X(1) ,Y(1) ,Z(1) 

COMMON/MAIN 1/NDIM 
N1M1=N1-1 
NN3=N3*NDIM 
DO  1 I=1,NN3,NDIM 
II=I+N1M1 
DO  1 J=I,II 
1 Z(J)=0.0 

ENTRY  MMULS 
NN3=N3*NDIM 
KK=0 

DO  10  K=1,N2 
DO  8 1=1. N1 

C1=X(I+KK)  , , , ,, 

IF(CI.NE.O.O)  CALL  VADD1 ( NN3, Cl , Z( I) ,Y(K) ) 
8 CONTINUE 

10  KK=KK+NDIM 

RETURN 
END 
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C 

C 


4 

5 

10 

6 
7 


SUBROUTINE  DIAG2(N,A,B,C1 ,C2) 

A = C1»B  + C2»I 

A,B  ARE  N»N  MATRICES;  I IS  N»N  IDENTITY  MATRIX 
DIMENSION  A(1) , B(1) 

COMMON  /MAIN1/  NDIM,  NDIM1 

NN=N»NDIM 

NM1=N-1 

11=1 

IF  (Cl  .EQ.  1.0)  GO  TO  10 
DO  5 J=1,NN,NDIM 
K=J+NM1 
DO  4 I=J.K 
A(I)=C1»B(I) 

Aai)=A(llT+C2 

II=II+NDIM1 

RETURN 

DO  7 J=1,NN,NDIM 

K=J+NM1 

DO  6 I=J,K 

A(I)=B(lj 

A(II)=A(II)+C2 

II=II+NDIM1 

RETURN 

END 
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SUBROUTINE  IDENT( N ,A , Cl ) 
DIMENSION  A(1) 
COMMON/MAIN 1 /NDIM , NDIM1 
NN=N*NDIM 
11=1 

DO  1 1=1, N 
DO  2 J=I,NN,NDIM 
2 A(J)=0.0 

A(II)=C1 

1 II=II+NDIM1 

RETURN 
END 


8? 


oo 
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SUBROUTINE  EQUATE(A,B,NR,NC) 
A=B 

MATRIX  EQUATE 
DIMENSION  A(1),  B(1) 

CALL  MSCALE(A,B,NR,NC, 1 .0) 

RETURN 

END 
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SUBROUTINE  MSCALE(A,B,NR,NC,C1 ) 
C A=C1*B 

c a"and  b may  be  equivalent 

DIMENSION  A(1) , B(1) 

COMMON  /MAIN1/  NDIM 
NN=NC*NDIM 

IF  (Cl  .EQ.  1 .0)  GO  TO  10 

IF  (Cl  .EQ.  0.0)  GO  TO  20 

IF  (Cl  .EQ.  -1 .)  GO  TO  30 

DO  5 1=1, NR 
DO  5 J=I,NN,NDIM 
5 A(J)=C1»B(J) 

RETURN 

10  DO  15  1=1, NR 

DO  15  J=I,NN,NDIM 
15  A(J)=B(J) 

RETURN 

20  DO  25  1=1 , NR 

DO  25  J=I,NN,NDIM 
25  A(J)=0.0 

RETURN 

30  DO  35  1=1, NR 

DO  35  J=I,NN,NDIM 
35  A(J)=-B(J) 

RETURN 

END 
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FUNCTION  D0T(NR,A,B) 

DOUBLE  PRECISION  DDT1,  DBLE 
DIMENSION  A(1),B(1) 
DDTIrO.ODO 

IF  (NR  .LE.  0)  GO  TO  2 
DO  1 1=1, NR 

1 DDT1=DDT1+DBLE(A(I)*B(D) 

2 DOTzDDTI 
RETURN 
END 
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FUNCTION  D0T2(NN,A,B) 

DOUBLE  PRECISION  DDT2,  DBLE 
DIMENSION  A(1) ,B(1) 

COMMON  /MAIN1/  NDIM 
DDT2=0.0D0 

IF  (NN  .LE.  0)  GO  TO  2 
DO  1 I=1,NN,NDIM 

1 DDT2=DDT2+DBLE(A(I)*B(D) 

2 D0T2=DDT2 
RETURN 
END 
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FUNCTION  D0T3(N,A,B) 

DOUBLE  PRECISION  DDT3,  DBLE 
DIMENSION  A(1) ,B(1) 

COMMON  /MAIN1/  NDIM 
DDT3=0.0D0 

IF  (N  .LE.  0)  GO  TO  2 
11=1 

DO  1 1=1, N 

DDT3=DDT3+DBLE( A(II)*B(I) ) 

1 II=II+NDIM 

2 DOT3=DDT3 
RETURN 
END 
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SUBROUTINE  VADD(N,C1 ,A,B) 
DIMENSION  A(,1)  ,B(1) 

DO  1 1=1, N 

1 A(I)=A(I)+C1*B(I) 

RETURN 

END 
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SUBROUTINE  VADD1( 
DIMENSION  A(1),B( 
COMMON/MAIN 1/NDIM 
DO  1 Ir1.NN.NDIM 
1 A(I)=A(I)+d»B(I) 

RETURN 
END 


NN,C1 ,A,B) 
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1 


5 


8 


13 


SUBROUTINE  VSCALE(X,Y,N,C1 ) 
DIMENSION  X(1) ,Y(1) 

L=0 

IF(c1 .EQ.1 .0)  GO  TO  5 
IF  01 .EQ.O.Oj  00  TO  8 
IF(C1 .EQ.-1 . J GO  to  13 
L=L+1 

X(L)=C1»Y(L) 

IF(L.LT.N)  GO  TO  1 

RETURN 

L=L+1 

X(L)=Y(L) 

IF(L.LT.N)  GO  TO  5 

RETURN 

L=L+1 

X(L)=0.0 

IF(L.LT.N)  GO  TO  8 

RETURN 

L=L+1 


X(L)=-Y(L) 
IF(L.LT.N)  GO  TO  13 
RETURN 
END 
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SUBROUTINE  VMAT1(A,X,N1 ,N2,Y) 
C Y=AX 

DIMENSION  A(1),X(1),Y(1) 

C0MM0N/MAIN1/NDIM 

DO  1 1=1, N1 

Y(I)=0.0 

II=I 

DO  1 J=1,N2 
Y(I)=Y(i5+A(II)*X(J) 

1 II=II+NDIM 

RETURN 
END 
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SUBROUTINE  VMAT2(Z,A,X,N1 ,N2,Y) 
C Y=Z+AX 

DIMENSION  A(1) ,X(1) ,Z(1) ,Y(1) 
C0MM0N/MAIN1/NDIM 
DO  1 1=1, N1 
Y(D  = Z(l) 

II=I 

DO  1 J=1,N2 
Y(I)=Y(I)+A(II)»X(J) 

1 II=II+NDIM 

RETURN 
END 
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FUNCTION  XGAIN(TH,XM,XS) 

DIMENSION  A(5) 

DATA  A/.2258368,-.2521287,1 .259695,-1 . 287822 9406461 / 
IF  (TH.GT.O.)  GO  TO  2 
XGAIN=1 .0 
RETURN 

2 Y = XM 

NS=2 

IF(XS.LT. 1 .0E-10)XS=1 .OE-10 

ifCy.eq.o.)  NS=1 
ANS=0. 

RMS=XS**2+XM**2 
DO  1 1=1, NS 
Z=.707*(TH+Y)/XS 
TEMP=EXP(-Z**2) 

X=1 ./(I .+.327591*ABS(Z)) 

P=X*(U(A(5)*X+A(4))»X+A(3))*X+A(2))*X+A(1))*1 .128379 
ERF=1 .-P*TEMP 
IF  (Z.LT.O.)  ERF=-ERF 

ANS=ANS+(RMS+TH*Y)*(1 .-ERF)-XS*Y*TEMP* .7975 
1 Y=-Y 

XGAIN=ANS/RMS/FLOAT(NS) 

IFCXGAIN.LT. 1 .E-6)  XGAIN=1.E-6 

RETURN 

END 
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SUBROUTINE  MATI0(X,NR,NC,I0) 

BATCH  ORIENTED  MATRIX  I/O 
10=1  INPUT  ONLY 

10=2  INPUT  AND  OUTPUT 

10=3  OUTPUT  ONLY 
10=4  PUNCH 

DIMENSION  X(1) 

COMMON  /MAIN1/  NDIM 

COMMON  /INOU/  KIN,  KOUT,  KPTR,  KPUNCH 

JEND=NC»NDIM 
GO  TO  (5,5,20,40)  10 

c»»»»»»»input 

5 DO  10  1=1, NR 

READ  (KIN, 1000)  (X(IJ),  U = I , JEND,NDIM) 

10  CONTINUE 

IF  (10  .EQ.  1)  RETURN 

C»»»»»»»0UTPUT 

20  DO  30  1=1, NR  , , , 

WRITE  (KOUT, 2000)  (X(IJ),  U=I , JEND,NDIM) 

30  CONTINUE 

RETURN 

C»»»»»»»puNCH 

40  DO  50  1=1. NR  . , , . 

WRITE  (KPUNCH, 3000)  (X( IJ) ,IJ=I , JEND,NDIM) 

50  CONTINUE 

RETURN 

1000  FORMAT  (8E10.0) 

2000  FORMAT  UH  ,1P10E13-3) 

3000  FORMAT  (1P8E10.3) 

END 
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SUBROUTINE  VECTIO(X, N, 10) 

BATCH  ORIENTED  VECTOR  I/O 
10=1.  INPUT  ONLY 
10=2  INPUT  AND  OUTPUT 
10=3  OUTPUT  ONLY 
10=4  PUNCH 

DIMENSION  X(1) 

COMMON  /INOU/  KIN,  KOUT,  KPTR,  KPUNCH 
GO  TO  (10,10,20,40)  10 

C*******iNPUT 

10  READ  (KIN, 1000)  (X(I),  1=1, N) 

IF  (10  .EQ.  1)  RETURN 

C*******0UTPUT 

20  WRITE  (KOUT, 2000)  (X(I),  1=1, N) 

RETURN 

c*******punch 

40  WRITE  (KPUNCH, 3000)  (X(I),  1=1, N) 

RETURN 

1000  FORMAT  (8E10.0) 

2000  FORMAT  (1H  .1P10E13.3) 

3000  FORMAT  (1P8E10.3) 

END 
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SUBROUTINE  PAGEFD(KFIL,KOUNT) 

C WRITES  KOUNT  FORMFEEDS  (1  IN  COL  1)  ON  FILE  KFIL 

ENTRY  FORMED 

100  IF  (KOUNT. LE.O)  RETURN 

DO  200  1=1, KOUNT 
WRITE  (KFIL, 1000) 

1000  FORMAT  (1H1) 

200  CONTINUE 


RETURN 

END 
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C 

10 

1000 


SUBROUTINE  DAYTIM(KFIL) 

WRITES  THE  DATE  AND  THE  TIME  ON  FILE  KFIL 

CALL  TIME(LTIME) 

CALL  DATE(LDATE) 

WRITE  (KFIL, 1000)  LDATE,  LTIME 
F0RMAT(1H  ,A10,2X,A10) 

RETURN 

END 
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